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This  thesis  is  intended  as  a source  document  for  further 


research  in  system  identification  processes.  Specifically, 
this  research  concentrated  on  efficiently  finding  the  inverse 
of  confluent  (i.e.,  repeated  eigenvalue)  Vandermonde  matrices, 
which  is  an  important  part  of  the  component  matrix  approach 
tc  system  identification.  However,  in  order  to  provide  a 

i 

broader  understanding  of  system  identification  in  general, 
and  the  component  matrix  approach  in  particular,  the  back- 
ground and  theory  sections  are  more  comprehensive  than  might 
otherwise  be  needed.  Several  supplementary  mathematical 
developments  and  some  example  problems  are  included  to  aid 
this  broad  understanding,  but  they  are  placed  in  appendices  to 
preserve  the  continuity  of  the  text,  and  referenced  as  appro- 
priate. A list  of  special  symbols  and  abbreviations  is  also 
included,  to  give  the  reader  a ready  reference  for  those  that 
may  be  unfamiliar. 
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Abstract 


Efficient  inversion  of  a 2n  x 2n  Vandermonde  matrix  is 
a key  requirement  of  a new  algorithm  developed  by  J.  Gary 
Reid  (Ref  10,  11)  for  parameter  identification  in  linear 
time-invariant  systems.  It  has  features  very  desirable  for 
application  to  large  systems  with  many  unknown  parameters, 
such  as  adaptive  flight  control  systems. 

A generalized  algorithm  for  inverting  the  Vandermonde 
matrix  was  proposed-  by  F.  G.  Cs&ki  (.Ref  1)  . It  was  chosen 
for  study  because  its  structure  parallels  that  of  Reid's 
algorithm.  Csaki's  algorithm  was  coded  into  a computer  sub- 
routine called  "VANINV,"  and  43  eigensystems  were  used  to 
test  its  computational  accuracy  and  efficiency.  Four  po- 
tential problem  areas  tested  were;  Cl)  large  system  orders, 
C2)  eigenvalues  near  each  other,  (.3)  eigenvalues  near  zero, 
and  (_4)  eigenvalues  with  large  magnitude  differences.  For 
a comparison,  the  Vandermonde  matrix  for  each  system  was 
also  inverted  using  routines  from  the  International  Mathe- 
matical & Statistical  Library  (IMSL) . 

Test  results  indicated  that  VANINV  is  not  quite  as  fast 
or  as  accurate  as  the  IMSL  routines.  However,  the  tests  were 
limited  to  systems  with  real  distinct  eigenvalues  because  the 
IMSL  routines  cannot  handle  other  types.  VANINV  in  its 
present  form  can  handle  systems  having  any  combination  of 
complex  and/or  repeated  eigenvalues.  Therefore,  recommen- 
dations for  further  research  and  several  possible  means  of 
improving  VANINV  are  outlined. 

ix 


INVESTIGATION  OF  INVERSE  VANDERMONDE  MATRIX 
CALCULATION  FOR  LINEAR  SYSTEM  APPLICATIONS 


Introduction 


Fundamental  quantities  needed  for  iterative  computation 
of  a quasi-linear  estimate  of  unknown  parameters  in  a linear 
dynamic  system  are  the  "parameter  sensitivity  variables." 
Finding  these  variables  is  a basic  part  of  a "system 
identification  process,"  as  is  shown  in  the  background  sub- 
section on  system  identification. 

Section  II  develops  the  process  for  finding  the  parameter 
sensitivity  variables  and  shows  that  a basic  part  of  the 
process  is  the  inversion  of  a 2n  x 2n  Vandermonde  matrix. 

The  development  of  an  efficient  computer  program  to  invert 
the  Vandermonde  matrix,  the  major  emphasis  of  this  thesis 
effort,  is  presented  in  Section  III.  Section  IV  explains 
some  test  procedures  used  for  determining  the  accuracy  and 
efficiency  of  the  program.  The  concluding  two  sections  give 
the  results  of  the  research  and  recommendations  for  further 
research. 

As  an  aid  in  understanding  the  concepts  and  processes 
presented  in  the  theory  sections  of  the  thesis,  several 
example  problems  and  supplementary  mathematical  develop- 
ments are  worked  out.  However,  because  of  their  length 
they  are  placed  in  appendices  and  referenced  as  appropriate. 
The  computer  program  developed  to  invert  the  Vandermonde 


matrix  and  the  programs  used  to  test  its  computational 
speed  and  accuracy  are  also  placed  in  appendices. 


would  suffice.  However,  control  systems  for  other  than 
very  simple  applications  are  not  linear  or  time-invariants 
depending  on  operating  conditions,^  the  control  laws  must 
change  to  maintain  adequate  system  control.  For  example, 
the  Wright  Brothers  didn't  have  to  concern  themselves  with 
large  changes  in  operating  conditions.  Their  first  powered 
flight,  120  feet  in  12  seconds  (Ref  2:Vol  29:557;  14 sVol  7: 
3881,  was  not  even  as  fast  as  a man  can  run,  and  didn't 
last  long  enough  for  the  environment  to  change  much. 

On  the  other  hand,  engineers  designing  flight  control 
systems  for  modem  high  performance  aircraft  do  need  to 
concern  themselves  with  large  changes  in  operating  condi- 
tions. As  aircraft  flight  envelopes  are  extended  over 
broader  operating  ranges  (for  example,  a high  performance 

1In  this  study,  "operating  conditions"  is  used  to 
encompass  both  system  configuration  and  environment. 

Each  has  definite  effects  on  system  response,  and  may 
also  affect  the  other.  For  convenience  however,  the  term 
"operating  conditions"  is  used  to  refer  to  both  when  it  is 
not  necessary  to  distinguish  between  them. 


aircraft  may  be  required  to  have  a Vertical  or  Short  Take- 
off and  Landing  (V/STOL)  capability) , more  and  alternate 
loops  must  be  added  to  the  flight  control  system  to  handle 
all  the  various  operating  conditions.  The  increased  re- 
quirement for  control  system  flexibility  often  results  in 
an  overall  system  that  is  nonlinear  and/or  time-varying, 
and  to  maintain  adequate  control,  some  means  of  adjusting 
the  control  system  parameters  is  necessary.  This  capacity 
to  change  the  modes  of  the  control  system  according  to  the 
aircraft  flight  operating  conditions  is  the  essence  of  an 
"adaptive"  flight  control  system.  Of  course  every  time 
the  operating  conditions  change,  the  control  system  vari- 
ables must  change  accordingly.  Thus  the  control  problem 
becomes  one  of  identifying  the  system  parameters  (the 
process  is  called  "system  identification")  in  order  to 
determine  the  proper  system  model  and  match  an  appropriate 
control  law  to  it. 

The  flight  control  system  design  engineer  could  try 
to  think  of  a collection  of  representative  types  of  oper- 
ating conditions  which  would,  hopefully,  cover  the  con- 
tinuum of  possible  operating  conditions.  A control  loop 
could  then  be  designed  to  handle  each  different  represen- 
tative operating  condition.  That  would  require  an  a priori 
identification  of  the  system  model  a nd  parameter  values 
and  a corresponding  collection  of  hardware  to  implement 
the  design.  Then  some  type  of  master  controller  would  be 
needed  to  choose  which  model  and  parameter  values  to  use 


for  the  given  operating  conditions. 

In  terms  of  an  entirely  analog  flight  control  system, 
the  increased  flexibility  required  to  handle  the  adaptive 
part  of  the  control  function  results  in  a confusing  mass 
of  hardware  and  interconnecting  links.  Also,  each  new 
piece  of  hardware  poses  a possible  reliability  problem. 
These  are  some  of  the  reasons  for  the  trend  toward  digital 
flight  control  systems.  Instead  of  having  a discrete 
piece  of  hardware  for  each  required  control  function,  a 
digital  controller  need  be  simply  reprogrammed  to  accom- 
modate the  new  control  system  parameters.  These  parameters 
provide  the  relation  between  the  system  inputs,  states, 
and  outputs  for  the  control  variable  calculations. 

The  overall  control  system  can  be  visualized  as  shown 
in  Figure  1.  Upon  starting  the  system,  some  nominal  system 
model  and  control  are  used.  However,  as  shown,  the  model 
must  change  to  account  for  significant  variations  in  the 
system  operating  conditions.  Referring  back  to  the  high 
performance  aircraft  mentioned  before,  an  example  of  a 
change  in  the  system  configuration  is  the  change  from  a 
vertical  takeoff  to  "conventional"  horizontal  flight.  An 
example  of  an  environmental  change  affecting  the  system 
model  is  the  change  in  air  density  with  altitude  (eg:  the 
air  at  40,000  feet  is  much  less  dense  than  at  sea  level, 
so  the  control  surface  actuator  gains  must  be  changed  ac- 
cordingly to  maintain  the  same  performance  relative  to 
pilot  inputs) . The  digital  control  system  can  be  given 


the  task  of  identifying  the  present  system  parameters, 
updating  the  system  model  with  them,  and  then  deter- 
mining what  adaptations  must  be  made  in  the  control  to 
achieve  the  desired  overall  performance.  The  first  two 
steps  are  what  distinguish  an  adaptive  control  system 
from  a fixed  control  system. 


System  Identification  Process  Backcrround 


A conceptually  simple  method  for  the  system  identifi 
cation  process  is  the  use  of  a first  order  Taylor  series 
approximation  to  update  the  system  model  according  to: 


where : 


y(t)  = physically  measured  system  output 


3y  Ct) 
aeT^T 


calculated  output  based  on  an  a priori 
model  Cor  in  an  iterative  process,  the 
last  updated  model) 

system  output  "sensitivity"  with  respect 
to  the  system  parameters 


A0  * change  in  the  system  parameters  from  those 
of  the  a priori  model  (or  the  change  since 
the  last  model  update,  for  am  iterative 
process) 

Notice  that  since  Ay  ■ A0,  Eq  (2)  is  equivalent  to 

y(t)  - y (t)  + Ay (t)  (3) 


which  appeals  to  the  intuitive  relation  "new  output  equals 
old  output  plus  the  change  in  the  output." 

It  is  seen,  then,  that  finding  these  sensitivity  variables 
is  a significant  part  of  the  system  identification  process. 
Therefore,  real-time  adaptive  control  applications,  such  as 
the  adaptive  flight  control  system  mentioned  earlier,  require 
an  efficient  means  of  computing  the  sensitivity  variables. 

However,  even  when  using  minimal  order  sensitivity  models, 
the  computational  requirements  can  be  quite  large  for  even 
modest  size  systems.  For  example,  the  complete  solution  of 
a sensitivity  system  of  order  n,  with  p unknown  parameters, 
may  require  the  solution  of  nCp+1)  coupled  linear  differential 
equations  CRef  8:1241.  The  basis  for  this  number  is  pre- 
sented in  Section  II. 

An  alternate  method  has  been  proposed  by  J.  G.  Reid 
in  his  dissertation  (Ref  11: Section  IIIl  which  may  have 
significant  computational  advantages  because  it  avoids  the 
need  to  solve  a collection  of  coupled  linear  differential 
equations  (CLDE} . This  method  is  explained  in  more  detail 
in  Section  II. 

However,  this  alternate  method  also  involves  the  in- 
version of  a 2n  k 2n  Vandermonde  matrix,  as  is  shown  in 
Section  II.  (Vandermonde  matrices  are  defined  in  Appendix  C, 
Eqs  (C2-2  through  C2-6)l . This  process  can  be  computationally 
quite  costly  (Ref  5:817).  Therefore,  a critical  part  of  the 
proposed  alternate  method  is  the  efficient  inversion  of  the 
Vandermonde  matrix.  Thus  the  effort  of  this  investigation 


concentrates  on  implementing  and  evaluating  a proposed 
method  (Ref  1:154-157}  of  efficiently  computing  the  in- 
verse Vandermonde  matrix.  The  inversion  method  used  forms 
polynomials  from  the  system  eigenvalues  and  uses  the  poly- 
nomials to  determine  the  elements  of  the  inverse  Vander- 
monde matrix.  Section  III  amplifies  this  method  and  shows 
its  parallel  computation  structure.  The  significance  of 
this  parallel  structure  is  explained  in  terms  of  its  poten- 
tial application  to  on-line  flight  controllers.  Section  IV 
explains  the  tests  used  to  evaluate  the  computation  time 
and  accuracy  of  the  computer  subroutine  written  to  implement 
the  method.  The  subroutine  was  called  VANINV.  Sections  V 
and  VI  give  the  test  results  and  resulting  conclusions  and 
resulting  recommendations  for  further  research. 


II  Parameter  Sensitivity  Calculations 


Section  I showed  that  calculation  of  the  parameter 
sensitivities  is  an  essential  part  of  the  system  identifi- 
cation process.  In  this  section,  a theoretical  comparison 
is  made  between  the  coupled  linear  differential  equations 
(CLDE)  approach  to  finding  parameter  sensitivities  and  the 
alternate  component  method  (CM)  approach  proposed  by  Reid. 
Since  the  CLDE  approach  is  already  developed  and  documented 
(eg:  Ref  13)  , it  is  merely  outlined  here  to  serve  as  a 
reference  base  for  the  CM  approach. 

Although  the  CLDE  method  is  more  straightforward,  it 
has  several  drawbacks  (Ref  8:124).  Perhaps  the  most 
obvious  difficulty  is  the  computational  burden.  As  is 
illustrated  later  in  this  section,  the  solution  of  the 
sensitivity  system  requires  solving  as  many  as  n(p+l) 
coupled  linear  differential  equations,  where  n is  the  sys- 
tem dimension  and  p is  the  number  of  parameters  whose 
sensitivities  are  to  be  found.1  Attempts  to  reduce  this 
computational  burden  with  low-order  sensitivity  models 
have  resulted  in  some  other  problems,  such  as:  (1)  require- 
ment for  special  forms  of  the  plant  matrix,  (2)  numerical 
sensitivity  of  calculations  performed  with  the  reduced 

1This  task  grows  increasingly  formidable  as  n and  p 
increase.  Each  new  variable  not  only  adds  another  CLDE, 
it  may  contribute  to  the  complexity  of  any  or  all  of  the 
other  equations. 
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order  models,  and  (3)  difficulty  of  obtaining  the  original 
state  variable  sensitivities  from  the  calculated  low-order 
sensitivities.  In  addition,  physical  insight  is  impaired 
when  using  the  transformed  coordinate  systems  resulting 
from  the  low-order  models. 

Among  the  advantages  of  the  CM  method  are:  Cl)  poten- 
tial reduction  in  computational  burden,  (2)  retention  of 
physical  system  insight  Cthe  modes  are  visible  throughout 
the  calculation  procedure) , and  essentially  "free"  avail- 
ability of  the  eigenvalue  sensitivities  themselves 
(Ref  11:79-80) . The  entire  set  of  parameter  sensitivities 
may  be  found  from,  at  most,  2nr  "quadrature  integrals" 

(where  n is  the  system  dimension  and  r is  the  number  of 
control  inputs),  plus  some  matrix  operations  (Ref  8:133). 

The  last  part  of  this  section  shows  that  the  CM  method 
requires  the  inversion  of  a 2n  x 2n  Vandermonde  matrix. 

Since  this  inversion  is  potentially  very  costly.  Section 
III  is  devoted  to  investigating  the  possible  method  that 
was  referenced  in  the  introduction  for  efficiently  computing 
the  inverse.  At  this  point,  some  system  definitions  are 
needed . 

Most  control  systems  other  than  very  simple  ones  are 
not  time- invariant.  However,  assuming  that  the  variations 
are  slow  w.r.t.  the  identification  process  (as  would  be  the 
case  for  flight  control  systems) , the  system  may  be  considered 
to  be  pseudo  time- invariant  (or  "quasi-static") . Thus,  for 
this  discussion,  the  linear  time-invariant  state  system  is 
used: 


] 


H 


J 

j 


10 


i 


( ) 


x (t ) 

= A ( 0) 

x (t) 

+ B (0)  u (t) 

(4a) 

y (t) 

= C ( 0 ) 

x (t) 

+ D(0)u(t) 

(4b) 

The  relation  of  the  initial  conditions  to  the  parameter 
sensitivities  is  defined  by  the  vector 


x(0)  = <j>  ( 0) 


(5) 


This  shows  that  the  system  matrices  and  the  initial  states 
are  a function  of  the  "parameter  vector,"  0,  of  which 
there  are  p elements,  designated  0^  i = 1,  2,  ...  p.  The 
nominal  value,  0Q,  of  the  parameter  vector  is  used  for  all 
evaluations.  Furthermore,  it  is  assumed  that  A,  B,  C,  D 
and  <j>  are  real,  bounded  and  continuously  differentiable 

w.r.t.  the  parameter  component  0^  at  its  nominal  value. 

The  "state  sensitivities"  of  this  system  are  defined 

as  the  partial  derivative  of  the  system  state  fector  w.r.t. 
the  parameter  components  0^,  and  are  written  as 


3x  ( t ; 0) 


(6) 


Coupled  Linear  Differential  Equations  Method 

With  the  assumptions  made  in  the  introduction  for 
this  section,  the  system  parameter  sensitivities  may  be 
found  by  augmenting  the  original  system  with  the  desired 
state  sensitivities  (from  Eq  (6) ) as  follows  (Ref  8:124}: 
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where  q is  the  system  output  dimension  and  r is  the  system 
input  (control)  dimension.  This  augmentation  gives  a 
"sensitivity  system"  which  contains  all  the  original  states 
plus  the  desired  parameter  sensitivities  expressed  as 
states.  From  Eq  (7)  it  is  evident  that  the  complete 
solution  of  this  sensitivity  system  involves  solving  up 
to  n(p+l)  coupled  linear  differential  equations. 
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Component  Matrix  Method 


The  CM  method  of  finding  the  parameter  sensitivities 


seems  at  first  to  be  much  more  indirect  than  the  CLDE  method. 


However,  after  all  the  development  of  the  theory  is  complete f 
it  turns  out  that  some  steps  can  be  eliminated  because  they 


serve  only  as  convenient  dividing  points  for  the  problem 


when  solved  manually.  In  a computer  solution  process,  once 


the  inverse  Vandermonde  matrix  is  found,  all  the  sensi- 


tivities of  the  parameters  and  the  original  system  can  be 


determined  using  the  elements  of  the  matrix. 


Eq  (4b)  , the  system  output  equation,  has  a time  domain 


solution  (Ref  12  :373) 


y(t)-CeACt"to)xCt0)+y'eA(t“TlBu  (tldx+D  u (t)  CIO) 


For  simplicity  of  discussion,  the  initial  time,  tQ,  and 


the  feed  forward  matrix,  D,  may  be  set  to  zero.  Since  the 


system  is  time- invariant , the  output  matrix,  C,  can  be 


taken  outside  of  the  integral.  Thus  Eq  (1)  becomes 


t 

y (t)  * CeAtx(o)  + cy  eA(t-T)Bu  (T)dx  (11) 


The  state  transition  matrix  (STM) , e , is  not  a function 


of  the  variable  of  integration,  so  it  may  be  removed  from 


the  integral.  Therefore 


y(t)  « CeAtx(o)  + ceAty*< 


At  / e Atb  u (t) dt 


Notice  that  the  integral  is  now  only  a function  of  one 


variable  and  is  therefore  simple  to  compute.  Such  single 


variable  integrals  are  called  "quadrature  integrals." 

Notice  that  there  are  a total  of  nr  quadrature  integrals 

-It  , 

since  e is  an  n x n matrix  and  B is  an  n x r matrix. 

Eq  (2)  established  the  . eed  to  find  the  system  out- 
put parameter  sensitivities,  so  differentiating  Eq  (12)  w.r.t. 
the  parameter  vector,  8,  yields 
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(13) 


Note  that  none  of  the  parameter  sensitivity  deriva- 
tives in  Eq  (13)  are  functions  of  time  except  that  of  the 
state  transition  matrix: 


At  _ 3 


.At 


(i) 


38. 


(14) 


The  sensitivity  derivative  of  the  state  transition  matrix 
is  potentially  very  costly  to  obtain  because,  in  the  form 
shown  in  Eq  (14) , it  must  be  reevaluated  for  each  desired 
point  in  time. 

It  can  be  shown  (see  Appendix  C,  Section  1) , that  the 

sensitivity  derivative  of  the  original  state  transition 
At 

matrix,  • is  contained  in  the  evaluation  of  the  state 

transition  matrix  of  the  augmented  sensitivity  system,  i.e: 
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2n*2n 


By  an  application  of  the  Cayley-Hamilton  theorem  (see 
Appendix  C,  section  2)  the  state  transition  matrix  of  the 
augmented  sensitivity  system  is  found  to  be  a summation 
of  products  of  the  remainder  polynomial  coefficients  and 
powers  of  the  sensitivity  system  matrix: 

j-1 

a..(t) 
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that  is:  the  power  of  a matrix  can  be  brought  inside  and 
applied  to  each  element  (Ref  11:64;  3:241-248).  There- 
fore applying  Eqs  (15)  and  (17)  to  Eq  (16)  results  in 


ajl 


Eq  (18)  shows  that  the  scalar  multiplication  and  sum- 
mation operations  can  be  applied  separately  to  each  section 

15 


’ r„  ?•  - , : v?:  • 


At  this  point,  it  is  still  necessary  to  find  the  partial  deri- 
vative of  2n  powers  of  the  original  system  matrix,  A,  and 
solve  2n  differential  equations  to  find  the  remainder  poly- 
nomial coefficients,  a^  (t) . The  sensitivity  problem  ap- 
pears to  have  become  even  more  complicated. 

However,  Reid  shows  (Ref  11:70)  that  the  partial  deri- 
vative of  the  system  state  transition  matrix  can  be  found 
from  a summation  of  products  of  partial  derivatives  of  the 
system  "component  matrices"  and  the  system  modes: 


Eq  (20)  is  a simplification  of  Reid's  formulation  and 
applies  for  the  case  of  distinct  eigenvalues.  It  is  used 
here  for  clarity. 


Eq  (20)  can  be  extended  to  the  augmented  system: 


of  the  partitioned  state  transition  matrix  of  the  aug- 
mented sensitivity  system.  So  the  partial  derivative 
of  the  original  state  transition  matrix,  which  is  needed 
to  find  the  parameter  sensitivities  according  to  Eq  (13) , 


can  be  evaluated  simply  as 


(i)  39^ 


(t) 


It  is  evident  from  Eq  (21)  that  for  every  mode,  e^jfc,  in 
the  original  system,  there  is  an  additional  mode,  te*^, 
in  the  augmented  system.  However,  the  additional  modes 
occur  only  in  the  sensitivity  derivative  of  the  STM.  Thus 
the  STM  of  the  original  matrix  is  preserved  by  the  CM 
method.  Only  the  original  modes  appear  in  the  STM  for 
the  original  system,  e.g.: 

eAt  ’ SUj,o)eXjt  <22) 

3-1 

i 

Thus  the  CM  method  is  seen  to  have  a physical  insight  ad- 
vantage in  comparison  to  the  CLDE  method:  the  modes  of 
the  system  are  clearly  visible  throughout  the  CM  solution 
process,  whereas  usually  they  are  not  in  the  CLDE  process. 
But  there  is  now  a total  of  2n  modes,  so  the  Vandermonde 
matrix  which  must  be  inverted  to  find  the  component 
matrices  for  the  sensitivity  system  is  2n  * 2n.  (See 
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Appendix  C,  Section  3 for  a discussion  of  how  the  com- 
ponent matrices  are  formed.)  And,  a total  of  2nr  quad- 
rature integrals  must  be  found. 

This  section  theoretically  compared  the  CLDE  method 
of  computing  system  parameter  sensitivities  to  the  al- 
ternate CM  method  proposed  by  Reid.  It  showed  that  when 
the  state  system  is  augmented  to  include  the  parameter 
sensitivities,  the  CLDE  method  requires  the  solution  of 
up  to  n(p+l)  coupled  linear  differential  equations.  To 
reduce  this  computational  burden,  reduced  order  models 
are  sometimes  used,  but  the  required  coordinate  transfor- 
mations often  impair  physical  insight  into  the  system  and 
make  it  difficult  to  obtain  the  original  system  sensi- 
tivities from  the  reduced  order  sensitivities.  The  alter- 
nate procedure,  the  CM  method,  was  shown  to  retain  the 
physical  significance  of  the  system  variables  throughout 
the  calculation  process,  and  to  exchange  the  formidable 
task  of  solving  n(p+l)  coupled  linear  differential 
equations  for  2nr  much  simpler  quadrature  integrals,  plus 
some  matrix  operations.  However,  the  potentially  costly 
process  of  inverting  a 2n  x 2n  Vandermonde  matrix  was 
also  shown  to  be  part  of  the  CM  method.  Thus  Section  III 
deals  with  the  subject  of  efficiently  inverting  the  Vander- 
monde matrix. 
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Ill  Vandermonde  Matrix  Inversion 
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Preceding  sections  established  the  need  to  calculate 
parameter  sensitivities  as  part  of  a system  identification 
process.  Two  methods  for  finding  those  parameter  sensi- 
tivities were  compared:  the  CLDE  method  and  the  CM  method. 

The  latter  was  chosen  for  further  study  based  on  theoretical 
comparisons  in  Section  II.  However,  a key  issue  associated 
with  this  method  is  the  need  to  invert  a 2n  x 2n  Vandermonde 
matrix.  This  section  addresses  a number  of  problems 
associated  with  inverting  the  Vandermonde  matrix,  and 
investigates  a generalized  method  of  Vandermonde  inversion 
proposed  by  F.  G.  Cs£ki  (Ref  1:154-157).  A special  feature 
of  Cs&ki's  generalized  algorithm  is  also  discussed.  That 
feature  is  that  the  algorithm  can  be  implemented  on  a 
parallel  architecture  processor,  which  for  a real-time 
flight  control  system  in  particular,  is  significant. 

Potential  Problem  Sources 

There  are  basically  five  sources  of  problems  asso- 
ciated with  computing  the  inverse  of  the  Vandermonde  matrix, 
all  of  them  related  to  the  characteristics  of  the  eigenvalues. 
They  are: 

(1)  Complex  eigenvalues 

(2)  Repeated  eigenvalues 

(3)  Eigenvalues  close  to  each  other 
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(4) .  Eigenvalues  close  to  the  origin 

(5)  Combinations  of  very  large  and  very  small 
eigenvalues 

Each  of  these  potential  problems  is  addressed  in  turn. 

Complex  Eigenvalues.  Vandermonde  matrices  formed 
from  complex  variables  may  be  handled  three  ways:  (1) 
conversion  of  the  complex  matrix  to  a real  one  (Ref  9: 
Problem  6.8),  (2)  use  of  two  separate  sets  of  real  vari- 
ables for  all  calculations — one  set  for  the  real  parts 
and  the  other  for  the  imaginary  parts  of  the  complex 
variables,  and  (.3)  use  of  complex  variables  directly. 
Investigation  of  Cs£ki's  generalized  algorithm  was  de- 
sirable because  of  its  parallel  implementation  capability, 
as  is  discussed  later  in  this  section.  The  algorithm  has 
an  inherent  capability  for  handling  complex  eigenvalues, 
so  conversion  of  the  complex  Vandermonde  to  a real  one  was 
unnecessary.  Few  control  systems  have  entirely  complex 
eigenvalues,  so  with  this  in  mind  the  second  method  of 
handling  complex  variables  was  chosen.  In  this  way,  tests 
for  complexity  can  be  performed,  and  for  calculations 
involving  only  real  quantities,  the  complex  part  of  the 
calculations  can  be  omitted  to  save  computation  time. 

Repeated  Eigenvalues.  For  systems  with  repeated  eigen- 
values, the  derivative  equations  (Eqs  (C2-4,  C2-5) ) must  be 
used  in  the  formation  of  the  Vandermonde  matrix.  If  the 
Vandermonde  matrix  was  not  defined  this  way,  any  system 
having  repeated  eigenvalue  would  have  a singular,  and  thus 
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uninvertable,  Vandermonde  matrix.  CsAki's  generalized  al- 

c 

gorithm  never  forms  the  Vandermonde  matrix,  but  incorporates 
the  derivative  feature  directly  in  the  formation  of  the 
inverse  Vandermonde  matrix. 

Nearly  Identical  or  Nearly  Zero  Eigenvalues.  Eigen- 
values that  are  close  to  each  other  or  close  to  the  origin 
may  cause  numerical  instability  because  of  the  finite  word 
length  of  digital  computers.  Algorithms  using  differences 
of  eigenvalues,  which  Ca£ki's  generalized  algorithm  does, 
are  especially  susceptible  to  this  problem.  The  accuracy 
tests  of  this  thesis  therefore  include  some  systems  with 
close  eigenvalues  and  some  nearly  zero  eigenvalues. 

Eigenvalues  With  Large  Magnitude  Differences.  Combi- 
nations of  eigenvalues  with  large  differences  in  magnitude 
may  cause  scaling  problems  in  the  calculations.  Some 
tests  to  investigate  this  problem  are  also  included  in 
the  thesis. 

Cs&ki's  Generalized  Algorithm 

F.  G.  Cs&ki  proposed  (Ref  1:154-157)  an  algorithm  for 
finding  the  inverse  of  a Vandermonde  matrix  directly  from 
a collection  of  eigenvalues  and  their  multiplicities.  The 
Vandermonde  matrix  itself  is  never  formed  by  the  algorithm. 
Since  the  algorithm  finds  each  element  of  the  inverse  Vender 
monde  matrix  independently,  it  may  be  implemented  on  a 
parallel  architecture  computer  to  increase  the  real-time 
( computation  speed.  This  is  discussed  in  more  detail  in  the 
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next  part  of  this  section.  Many  inversion  techniques  can- 
not be  implemented  in  parallel  fashion  because  they  do 
simple  row  or  volumn  operations  rather  than  operations 
involving  individual  elements. 

Cs&ki's  generalized  algorithm  partitions  the  inverse 
Vandermonde  matrix  into  m blocks,  VT,  (i  = 1,  2,  . . . . m)  , 
with  one  block  for  each  distinct  eigenvalue.  Each  block 
is  dimensioned  k^  rows  by  n columns,  where  k^  is  the  eigen- 
value multiplicity  and  n is  the  system  order.  The  Vander- 
monde inverse  may  thus  be  represented  by 


where 


W 


(23) 


(24) 


in  which  WjTti)  ( j - 1,  2, 
of  the  block 


k^xn 


. , k^)  are  the  row  vectors 
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With  this  formulation  and  nomenclature 


in  which  represents  the  element  in  the  j ' th  row  and 

3 *- 

i'th  column  of  the  i'th  block,  W..  The  fraction  terms  are 


D(s)  is  the  system  characteristic  equation,  which  may  be 
represented  as 


Eq  (.26). , the  i'th  "denominator  polynomial,"  is  therefore 
the  characteristic  equation  with  the  i'th  eigenvalue 
divided  out  completely.  The  terms  defined  by  Eq  (27)  are 
called  "truncated  polynomials"  and  are  formed  by  itera- 
tively reducing  by  one  the  power  of  each  term  of  the 
characteristic  equation,  Eq  (28) , and  deleting  the  right- 


most  term  until  a total  of  n truncated  polynomials  are  found 
For  example,  suppose  a third  order  system  has  the  charac- 
teristic polynomial 


The  n = 3 truncated  polynomials  are 


C30c) 


Uote  that  the  final  truncated  polynomial  is  always  the 
integer  one  because  of  the  form  of  D(s)  used.  Also 
notice  from  Eq  (25)  that  the  derivative  procedure  which 


for  the  elements  of  the  inverse  matrix 


In  order  to  draw  the  foregoing  together,  suppose  a 
sixth  order  system  has  three  eigenvalues,  having  multi- 
plicities of  three,  two,  and  one,  respectively.  Substi- 
tuting Eq  (25)  into  Eq  (24)  , and  the  result  into  Eq  (23) 
yields  the  following  structure  for  the  inverse  Vandermonde 


v"1  = w = 


A simple  numerical  example  of  the  application  of  this 
algorithm  is  in  Appendix  D,  Example  3. 

Parallel  Architecture  Processors  and  Cs£ki*s 
Generalized  Algorithm 

If  a large  problem  c an  be  broken  up  into  a number  of 
smaller  independent  ones,  parallel  processing  can  be  used 
to  reduce  the  real  process  time.  No  computational  ad- 
vantage is  realized,  but  since  several  calculations  can 
be  done  simultaneously,  the  clock  time  from  start  to 
finish  of  the  large  problem  is  less  than  if  each  calcu- 
lation were  done  serially.  This  is  particularly  impor- 
tant for  flight  control  applications  where  real  processing 
time  may  be  critical. 

For  a comparison  between  serial  and  parallel  process 
times,  suppose  a third  order  system  has  two  eigenvalues, 
one  with  a multiplicity  of  two.  A serial  processor  would 
perform  the  Vandermonde  inversion  with  Csiki's  algorithm 
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as  shown  in  Figure  2a.  If  each  calculation  takes  p seconds 
of  processing  time,  the  minimum  total  time  is  15p  because 
a calculation  cannot  begin  until  the  preceding  one  is 
finished.  Figure  2b  shows  the  same  calculations  as  a 
parallel  processor  could  perform  them.  Notice  that  the 
total  time  from  start  to  finish  of  the  group  of  calculations 
is  a third  of  the  amount  of  time  the  serial  processing  took 
because  the  parallel  processing  system  has  three  times  as 
many  processors  to  do  the  job.  More  information  on  parallel 
processing  is  found  in  Ref  4. 

Not  only  does  Csciki's  method  for  finding  the  inverse 
of  the  Vandermonde  matrix  lend  itself  to  parallel  processing, 
so  does  the  component  matrix  approach  itself,  as  Figure  3 
shows.  At  the  top  of  the  figure  is  9^  the  current  estimate 
of  the  system  parameter  vector.  From  that  estimate,  the 
system  plant  model,  a (A} , can  be  found,  and  from  it,  the 
current  eigenvalues,  ...,  An,  can  be  determined. 

As  shown,  each  eigenvalue  can  be  used  in  parallel  to  com- 
pute the  1 x 2nr  vector,  F,  of  quadrature  integrals,  while 
at  the  same  time  computing  the  elements  of  the  inverse 
Vandermonde  matrix  which  correspond  to  that  eigenvalue. 

At  this  point,  the  value  of  the  parallel  computation  struc- 
ture of  Csiki's  algorithm  is  readily  visible.  If  the 
algorithm  to  find  the  inverse  Vandermonde  matrix  did  not 
have  a parallel  structure,  none  of  the  calculations  using 
the  elements  of  the  inverse  Vandermonde  matrix  could  pro- 
ceed until  the  entire  matrix  was  computed  in  serial  form. 
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Parallel  Processing  Steps 
with  n Processors 


Figure  3.  Parallel  computation  structure  of  the  component  matrix  method 
of  system  identification.  (Adapted  from  slides  used  to 
present  Ref  10.) 


That  could  stretch  out  the  overall  computation  cycle  time 
considerably,  as  Figure  2 showed.  However,  with  the 
parallel  structure  of  CsSki's  algorithm,  each  eigenvalue 
branch  of  Figure  2 may  find  the  elements  of  the  inverse 
Vandermonde  matrix  pertaining  to  that  eigenvalue  and 
continue  with  the  next  calculations  immediately.  From 
the  calculated  inverse  Vandermonde  matrix  elements,  the 
vector  of  component  matrices,  Y,  and  matrix  of  component 
sensitivities,  G,  can  be  computed.  A change  in  the  system 
parameters,  A0,  can  then  be  computed  from  (Ref  10) 

Y =[F]*  [G]TA0  (.32) 

where  Y,  [F] , and  [G]  are  defined  as  in  Figure  3. 

It  is  now  possible  to  use  A0_  to  update  the  estimate  of  the 

A 

parameter  vector,  Q_.  This  closes  the  computational  loop, 
making  the  system  identification  process  iterative,  and 
thereby  capable  of  applications  in  on-line  flight  controllers. 

All  of  this  sounds  nice  in  theory.  The  immediate 
question  that  arises  concerns  whether  or  not  the  method 
actually  works  when  implemented  on  a computer  with  a prac- 
tical problem  like  limited  word  length  which  causes 
truncation  and  round  off  errors  in  the  calculations  per- 
formed. Therefore  Cs4ki's  algorithm  was  encoded  in  Fortram 
and  some  tests  were  run  to  evaluate  its  accuracy  and 
efficiency.  Those  tests  are  the  subject  of  the  next 
section. 
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Tests  Run  on  the  Coded  Version  of 


Csaki's  Algorithm 


Previous  sections  of  this  thesis  dealt  with  the  theory 
behind  a system  identification  process.  The  CM  method  for 
system  identifications  was  shown  to  have  theoretical  compu- 
tational advantages  over  the  more  straightforward  CLDE 
method.  However,  part  of  the  CM  method  is  the  inversion 
of  a 2n  x 2n  Vandermonde  matrix,  which  can  be  computa- 
tionally costly.  Csaki's  generalized  algorithm  for  finding 
the  inverse  Vandermonde  matrix  was  chosen  for  its  parallel 
computation  structure.  This  structure  gives  the  algorithm 
the  capability  of  operating  much  faster  in  real  time  than 
algorithms  which  must  perform  all  the  calculations  in  serial 
form.  Also,  the  particular  parallel  structure  of  Csaki's 
algorithm  is  very  compatible  with  the  parallel  structure 
of  the  CM  method  of  system  identification. 

Csaki's  generalized  algorithm  was  coded  in  Fortran  in 
a subroutine  called  "VANINV",  and  is  hereafter  referred  to 
as  such.  Appendix  A contains  descriptions  and  flowcharts 
of  VANINV  and  the  subroutines  associated  with  it.  Part  I 
of  Appendix  B contains  source  listings  of  all  the  sub- 
routines in  Appendix. A.  This  section  describes  tests  done 
subroutines.  Source  listings  of  the  test  programs,  though 
not  flowcharted,  are  found  in  Appendix  B,  Part  II. 
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Types  of  Tests 

Two  types  of  tests  were  performed  on  VANINV.  The  first 
was  actually  a combination  of  four  separate  tests  to  deter- 
mine VANINV' s computational  accuracy.  The  second  type  of 
test  measured  how  long  it  took  VANINV  to  find  the  total 
inverse  Vandermonde  matrix  for  a given  system.  These  two 
types  of  tests  were  run  on  a number  of  different  systems 
to  determine  what  happens  to  the  computation  time  and 
accuracy  when:  (1)  the  system  order  increases,  (2)  two 
eigenvalues  become  very  close,  causing  the  Vandermonde  matrix 
to  become  more  nearly  singular  and  therefore  harder  to  in- 
vert, C3)  one  eigenvalue  approaches  zero,  which  also  makes 
the  Vandermonde  matrix  more  nearly  singular,  and  (4)  the 
difference  between  the  magnitudes  of  the  largest  and  smallest 

eigenvalue  increases.  The  ratio,  C = | -x. | , in  (4)  is 

Amin 

called  the  "condition  index." 

In  order  to  obtain  a relative  measure  of  VANINV' s 
performance,  the  same  tests  with  the  same  systems  were  run 
using  a standard  matrix  inversion  routine  contained  in  the 
International  Mathematical  & Statistical  Library  (IMSL) . 

Early  tests  were  done  with  IMSL  subroutine  LINV1F,  but  it 
did  not  have  sufficient  accuracy,  so  the  final  tests  were 
done  using  IMSL  subroutine  LINV2F,  which  has  a capability 
for  iterative  improvement  of  the  solution. 

Accuracy  Tests 

The  first  accuracy  test  is  actually  only  a test  of  the 
IMSL  routine.  It  simply  involves  feeding  the  inverse 


31 


* 


Vandermonde  matrix  hack  into  the  routine,  and  comparing 
the  reinverted  matrix  to  the  original  Vandermonde  matrix. 

It  is  unlikely  that  the  accuracy  of  the  reinverted  matrix 
is  better  than  that  of  the  first  inversion.  Thus  a lower 
accuracy  bound  on  the  inversion  is  the  minimum  number  of 
decimal  places  in  the  reinverted  matrix  that  are  identical 
to  the  original  matrix. 

Originally  a second  accuracy  test  was  set  up  to 
multiply  the  Vandermonde  matrix  by  the  calculated  inverse 
and  subtract  the  identity  matrix: 

V v"1  - I = E (33) 

in  which  the  subscript  "c"  means  "calculated."  Theo- 
retically, E should  always  be  a null  matrix.  The  extent 
to  which  it  is  not  zero  gives  an  indication  of  the  ac- 
curacy of  the  calculation  of  the  inverse,  V~^ . However, 

E is  typically  composed  of  numbers  so  small  that  the  errors 
caused  by  the  computational  characteristics  of  the  algo- 
rithm cannot  be  distinguished  from  the  truncation  and 
round  off  errors  resulting  from  the  finite  word  length  of 
the  computer. 

The  following  analysis  may  be  made  of  the  problem 
(Ref  6:61-64).  What  is  really  desired  is  the  error 
between  the  true  inverse  and  the  actual  inverse,  defined 
as: 
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Ac  - A = Error 


A(A~-A-1)  = Aa“1-AA"1=Aa“1-I  = E (35) 

c c c 

from  Eq  (33).  Then  Eq  (34)  can  be  expressed 

(A~1_A_1)  = Error  = A-1(AA~1-I)  = A_1E  (36) 
c c 

Tests  with  matrices  having  known  inverses  show  that  the  abso- 
lute value  of  the  smallest  negative  exponent  in  the  elements 

of  A ^E  approximates  the  minimum  number  of  decimal  places 
in  A-1  that  are  accurate.  For  example,  if  the  elements  in 
A-1E  are  .xxxxxE-10,  the  elements  of  A”1  may  be  considered 
accurate  to  at  least  ten  decimal  places.  This  may  be  ex- 
pressed as  follows: 


a^j  be  the  elements  of  A~ 
e^  be  the  elements  of  E 

be  the  elements  of  A-1E 


Therefore 


k-1 


ik  ekj 


and  a bound  on  the  error  of  is  given  by 


6ij  i * 


Ii-i 


However,  aik  e A-  is  not  available;  so  use  yik  e a"  and 
use  the  following  approximation  for  the  error  bound: 


->9*'*' 


(39) 


Kji  ^iki  i*kj 


hoping  that  the  actual  error  of  the  element  is  bounded  by 

p . . : 

I “ik  ‘ Yikl  < I °.J  (40) 


When  this  test  was  applied  to  systems  with  known  solutions, 
Eq  (40)  predicted  inversion  accuracy  to  within  three 
decimal  places  for  systems  whose  condition  index  was  less 
than  about  100,  unless  there  were  two  eigenvalues  closer 
them  .2  units  to  each  other. 

The  third  test  used  was  simply  a summation  of  the 
elements,  £5.^,  of  the  check  matrix,  A^E,  from  test  two. 

It  was  an  attempt  to  quantify  in  a single  number  the  ac- 
curacy of  the  inversion  method  being  tested.  In  this  thesis 
it  is  referred  to  as  an  "Accuracy  Merit  Figure"  CAMF) . The 
problem  with  it  is  that  it  is  very  dependent  on  how  many 

elements  of  A 1 are  nonzero.  It  tends,  therefore,  to  be  a 
c 

cumulative,  rather  than  an  absolute  (i.e.:  "decimal  place") 
accuracy  test.  Note,  however,  that  the  AMF  can  be  used  in 
conjunction  with  the  largest  element  of  the  check  matrix  to 
carry  the  information  contained  in  the  matrix  in  two  numbers 
rather  than  the  n elements  of  the  matrix.  Tests  showed 
that  if  the  AMF  is  much  larger  than  the  largest  element  of 
the  check  matrix  (i.e.:  AMF  >>  (B^ ) max)  * it  indicates  that 
many  elements  of  the  matrix  are  nearly  the  same  size  as 
($ij)max  30(5  probably  all  the  elements  of  the  inverse  matrix, 
V-1,  are  accurate  to  the  same  number  of  decimal  places. 


However,  if  AMF  = CS^niax'  one  element  of  the  check  matrix 
is  dominant,  and  most  of  the  other  elements  of  V-1  are  more 
accurate  than  indicates.  In  this  situation, 

^ j IUOA 

(Bij)max  Is  definitely  a worst  case  error  indicator. 

Fourth  and  final  of  the  accuracy  tests  used  was  a 
subjective  comparison  of  the  inverse  Vandermonde  matrix 
calculated  by  VANINV  and  the  IMSL  routine.  A number  of 
the  systems  used  in  the  tests  had  known  solutions,  or  had 
repeating  decimals  in  the  solution.  The  decimal  place 
accuracy  of  each  inversion  was  therefore  estimated  on  the 
basis  of  these  characteristics  of  the  solutions. 

A total  of  43  different  systems  were  designed  in 
order  to  test  VANINV  and  the  IMSL  routines.  The  tests  were 
run  on  a CDC  6600  computer  using  single  precision  arithme- 
tic. Results  of  the  tests  are  given  in  the  next  section. 
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V Results  of  Tests 

The  total  of  43  test  systems  (see  Table  I)  which  were 
run  were  divided  among  the  following  five  application 
categories: 

1.  Computation  time  as  a function  of  system  order 
(Tests  1 through  43)  . 

2.  Accuracy  as  a function  of  system  order  (Tests  1 
through  11) . 

3.  Accuracy  as  a function  of  the  distance  between 
the  closest  two  eigenvalues  of  the  system, 
lXx”Xylmin  CTests  12  through  20). 

4.  Accuracy  as  a function  of  the  condition  index, 

Unaxl 

C = ■rvin  -|  CTests  21  through  27)  . 

* min ' 

5.  Accuracy  as  a function  of  increasing  condition 
index,  C,  as  the  smallest  eigenvalue,  l^m^nl» 
of  the  system  approaches  zero  CTests  28  through 
43)  . 

Table  II  contains  all  the  data  from  the  tests,  and 
the  results  are  explained  according  to  the  categories  above. 

Computation  Time 

Tests  1 through  11  indicate  that  the  computation  time 

for  finding  the  inverse  Vandermonde  matrix  is  proportional 

2 3 

to  a number  between  n and  n . This  is  to  be  expected 

2 

since  the  number  of  elements  to  be  computed  is  n and  some 
elements  take  more  calculations  than  others.  Figure  4 
shows  this  relation  (computation  time  between  n * 8 and 
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Table  I 
Teat  Systems 


Tast  No. 


System  Order 


Xj-  2.  +J0. 
X2— 2.  +J0. 


Xj-  2.5+JO. 
Xj— 2.5+JO. 

Xj-  2.  +J0. 


Xj-  3.  +J0. 
X j"“ 3 . +J0. 
X3-  2.5+JO. 


Eigenvalues 


3 - 1.  +J0. 

4 S-I.  +J0. 


— 2.  +J0. 

- l.S+JO. 
—1.5+JO. 


- 2.5+JO. 

- 2.  +J0. 
—2.  +j0. 


X#  - 1.  +J0. 
X,  —1.  +J0. 


X?  - 1.5+JO. 
x8  — 1.5+JO. 
X-  “ 1.  +J0. 


X,  - -. 5+JO . 


3-  3. 5+JO. 

Xj  - 2 .5+ jO . 

Xg  - 1.5+JO. 

Xu-  .5+JO. 

2— 3. 5+JO. 

X#  —2.5+JO. 

X^g— 1.5+JO. 

X3g»  -.5+JO. 

3»  3.  +J0. 

X?  - 2.  +J0. 

x31»  1.  +J0. 

4— 3.  +JQ. 

X#  »-2 . +J0. 

Xl2— !•  *30. 

3-  4.  +J0. 

Xj  - 3.  +J0. 

X,  - 2.  +J0. 

xl3-  1.  +J0. 

2— 4.  +J0. 

X6  —3.  +J0. 

X10*-2,  *J0* 

xl4— 1.  +J0. 

3-  3.S+JQ . 

X?  - 2.5+jO. 

X13-  1.5+JO. 

X15-  ,5+JO. 

4*-J.5*30. 

Xg  —2.5+JO. 

X12— 1.5+JO. 

Xls-  -.5+JO. 

3-  4 .5+JO. 

Xg  —3. 5+JO. 

Xn-  2.  +J0 . 

X3g—1.  +J0. 

l2— 4.5+JQ. 

X,  - 3.  +J0. 

X12— 2.  +J0. 

Xn-  .5+JO. 

Ij"  4.  +J0. 

Xg  —3.  +JO. 

X13-  1.5+JO. 

X^-  -.5+JO. 

4— 4.  +J0. 

x#  - 2.5+30. 

X14—1.5+J0. 

kj-  3. 5+ JO. 

X10-  2.5+JO. 

Xu"  1.  +J0. 

L-  *.  +J0. 

Xfi  - 4.  +J0. 

X1X«  2.5+JO. 

Xig-l.S+JO. 

2— 5.  +J0. 

X?  - 3. 5+JO. 

X12«-2.5+J0. 

X17-  1.  +J0. 

3-  4.5+JO. 

Xg  —3. 5+JO. 

X13*  2‘  *i0' 

X3g»-1.  +JO • 

g— 4.5+JO. 

X,  - 3.  +J0. 

X^g— 2.  +J0. 

X19“ 

j-  4.  +J0. 

X1Q— 3.  +J0. 

X15«  1.5+JO. 

X16«  -.5+JO. 

Table  I — continued 
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Table  II 
Test  Results 


Table  II — continued 


njHDnBM9ff£J7BDI 

■BaP^pr^Eiai 

■arcarcsBB 


fmm 

marca 


psaHSBi 

RUIN 

MJbRW 


u It. (iiis-1!  | ».»t*to- 


u lj.i4«io*11 


U Ij.ajxto-1*  t.»««io-11 


1.44XlO*L> 

l.MsM*11 

l.7«xt«-lJ 

i.JOxlO'1* 

t.ta«wu 

4.47*t0'15 

*.u«io-10 

l.iJMlO"15 

t By  subjective  observation. 
* By  second  accuracy  test. 


* IN BL  had  a terminal  error. 
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COMPUTATION  TIME  (SEC) 

2.00  4.00  6.00  8. 00  1Q.00  12.00  14.00 


2 5 

n = 16  is  proportional  to  about  n * ) and  also  shows  that 
the  IMSL  routines  tended  to  run  only  slightly  faster  than 
VANINV.  The  remaining  tests  were  all  done  with  fourth  order 
systems,  and  computation  times  remained  nearly  constant  in 
every  case  for  VANINV. 

Accuracy  as  a Function  of  System  Order 

Figure  5 shows  graphically  what  the  results  of  Tests  1 
through  11  indicate — accuracy  is  essentially  constant  at  the 
maximum  machine  capability  (14  decimal  places  for  a CDC  6600 
machine)  over  the  range  of  system  orders  tested. 

The  second  accuracy  bound  test  indicates  that  although 
the  full  14  decimal  place  accuracy  is  present  up  through 
the  14th  order  system,  the  probability  of  this  being  so  is 
diminished  as  the  system  order  increases.  This  may  be  ex- 
plained as  follows.  The  number  of  calculations  done  to  find 

each  element  of  the  inverse  Vandermonde  matrix  increases 

2 3 

approximately  to  the  n to  n power  with  system  order,  as 
was  shown  in  the  preceding  result.  With  the  increase  in 
the  number  of  calculations,  truncation  and  round  off  errors 
from  the  calculations  may  accumulate  sufficiently  to  ap- 
pear to  propagate  into  the  more  significant  decimal  places. 
Since  the  second  accuracy  test  uses  absolute  values  of 
error  quantities  in  determining  a lower  found  on  the  ac- 
curacy, the  test  is  actually  a calculation  of  the  worst 
possible  accuracy,  not  the  most  probable  accuracy.  This 
hypothesis  would  be  supported  if  these  tests  were  run  in 
double  precision  arithmetic  and  the  accuracy  tracked  along 
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at  the  higher  level.  The  reason  for  the  sudden  drop  in  the 
minimum  bound  for  VANINV  beginning  with  the  18th  order 
system  is  not  intuitively  obvious,  but  it  is  interesting  to 
note  that  it  seems  to  be  indicative  of  a trend  as  the  sys- 
tem order  continues  to  increase. 


Accuracy  as  a Function  of  |AX-1  |m^n 


From  Figure  6 it  can  be  seen  that  the  accuracy  for  both 
VANINV  and  the  IMSL  routine  drops  off  as  the  minimum  distance 
between  two  eigenvalues  goes  to  zero.  Neither  routine  is 
significantly  better  than  the  other,  as  the  results  of 
Tests  12  through  20  show. 

Apparently  the  second  accuracy  bound  test  does  not  work 
very  well  for  systems  having  two  eigenvalues  close  to  each 
other.  The  closer  the  eigenvalues  become,  the  less  realistic 
the  accuracy  bound  is. 

It  is  interesting  to  note  that  although  the  accuracy 

of  the  IMSL  routine  decreases  as  | A —A  |min  decreases,  the 

x y 

resulting  inverse  is  stable,  in  each  test  the  inverse  and 
the  inverse  of  the  inverse  have  about  the  same  accuracy,  so 
the  solution  of  the  inverse  matrix  is  termed  "stable"  for 
this  discussion. 


Tests  21  through  27  hold  the  minimum  eigenvalue, 

IX  , I,  and  the  closest  distance,  |X  -X„l  between  two 
1 min 1 1 x y min  | ^ 

eigenvalues  constant  as  the  condition  index,  C ■ — 


Accuracy  as  a Function  of  Condition  Index 


eweisiesnmepaPMiifMHw 


'WHIP 


increased.  The  re inversion  of  the  inverse  by  the  IMSL 
routine  indicates  that  the  inversion  accuracy  decreases  as 
the  condition  index  increases.  The  actual  inversion  ac- 
curacy could  not  be  determined  by  inspection  of  the  data, 
though  VANINV  and  the  IMSL  routine  gave  identical  results 
up  through  a condition  index  of  twenty-five  thousand.  Be- 
yong  that  point  both  inversion  routines  became  unstable, 
as  indicated  by  the  large  fluctuations  in  the  elements  of 
the  check  matrix  and  the  Accuracy  Merit  Figure.  Plots  were 
not  made  of  this  test  series  because  of  the  small  amount  of 
plotable  data. 


Accuracy  as  a Function  of  Increasing  Condition  Index 


Caused  by  an  Eigenvalue  Approaching  Zero 


These  tests  were  run  in  two  series.  The  first  series, 
Tests  28  through  37,  allowed  one  eigenvalue  to  go  to  zero, 
thereby  causing  the  eigenvalue  condition  index,  C,  to  in- 
crease. The  minimum  distance  between  any  two  eigenvalues 
was  held  constant.  As  seen  in  Figure  7,  the  accuracy  of 
both  inversion  routines  dropped  as  the  condition  number 
increased. 

In  the  second  series  of  tests,  numbers  38  through  43, 
the  minimum  distance  between  any  two  eigenvalues  was  again 
held  constant,  but  at  a larger  value  in  order  to  determine 
if  it  has  any  effect  on  the  inversion  accuracy.  The  accuracy 
could  not  be  determined  by  inspecting  the  inverse  matrices 
found.  However,  based  on  the  accuracy  of  the  reinverted 


COMPUTATION  ACCURACY  (DECIMAL  PLACES) 

2.00  4.00  6.00  8.00  10.00  12.00  14.00  16.00 

. . . > ■ > • i 


.00  1.80  2.60  3.40  4.20  5.00 

LOGARITHM  OF  CONDITION  INOEX 


Figure  7.  Logarithm  of  condition 

index  Uog|AnuU[l/|*minl>. 
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matrix,  the  accuracy  and  stability  of  the  IMSL  routines 
dropped  as  the  condition  index  grew  in  each  series. 

From  the  results  of  both  series  of  tests,  a high  con- 
dition index  resulting  from  am  eigenvalue  near  zero  is  not 
as  detrimental  to  inversion  accuracy  or  stability  as  a high 
condition  index  caused  by  a large  eigenvalue.  In  addition, 
the  combination  of  close  eigenvalues  and  a high  condition 
index  has  more  of  a degrading  effect  on  accuracy  than  either 
of  the  two  by  themselves.  It  may  also  be  observed  that  when 
the  eigenvalue  actually  reaches  zero,  thereby  driving  the 
condition  index  to  infinity,  the  accuracy  rises  sharply  back 
to  the  level  realized  for  small  condition  indices.  This  is 
probably  because  the  computer  can  represent  zero  exactly, 
whereas  numbers  near  zero  can  only  be  approximated. 

All  of  the  findings  of  the  tests  are  summarized  in  the 
several  conclusions  of  the  next  section.  In  addition, 
recommendations  are  given  for  possible  improvements  to 
VANINV  and  for  its  use  in  further  research  regarding  system 
identification  processes. 


VI  Conclusions  and  Recommendations 


Several  conclusions  were  drawn  from  the  results 
reported  in  Section  V.  Based  on  those  conclusions  it  is 
recommended  that  subroutine  VANINV,  which  implements 
Cs£ki's  generalized  algorithm  for  finding  the  inverse  of 
a Vandermonde  matrix,  be  refined  for  computational  speed 
and  tested  in  its  intended  application,  that  of  finding 
component  matrices  for  a system  identification  process. 

The  conclusions  and  the  recommendations  are  explained 
further  in  the  remainder  of  this  section.  The  method  of 
presentation  is  a listing  of  the  conclusions/recommendations, 
followed  by  some  comments  regarding  them. 

Conclusions 

The  following  conclusions  were  made: 

1.  The  IMSL  routines  are  generally  slightly  faster, 
and  about  as  accurate  or  somewhat  more  accurate 
than  VANINV  is  in  its  present  form. 

2.  The  accuracy  of  both  routines  is  reduced  by  high 
condition  indices,  particularly  if  the  high  con- 
dition indices  are  caused  by  large  eigenvalues, 
as  opposed  to  being  caused  by  one  eigenvalue  near 
zero. 

In  regard  to  computation  speed,  it  should  be  noted 
that  VANINV  is  compiled  to  optimization  level  two,  whereas 
the  IMSL  routines  are  only  compiled  to  optimization  level 
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one.  This  means  that  IMSL  could  possibly  run  significantly 
fastf  r than  VANINV.  However,  these  tests  were  restricted 
to  real  distinct  eigensy stems.  VANINV  has  the  capability 
of  computing  the  inverse  of  Vandermonde  matrices  with  com- 
binations of  repeated  and/or  complex  eigenvalues.  A fair 
test  of  this  property  would  necessitate  special  conversions 
to  turn  the  Vandermonde  matrix  into  a purely  real  matrix 
that  the  IMSL  routines  could  handle.  This  conversion  time 
would  then  be  considered  part  of  the  IMSL  computation  time. 

It  would  be  interesting  to  determine  VANINV' s accuracy 
if  double  precision  arithmetic  was  used.  Of  course,  for 
any  intended  application  a trade-off  analysis  should  be  per- 
formed between  any  increase  in  accuracy  and  the  cost  in 
terms  of  additional  memory  requirements. 

Recommendations 

It  is  recommended  that  VANINV  be  developed  further  and 
implemented  in  a component  matrix  approach  for  system  iden- 
tification. In  this  way,  VANINV' s accuracy  can  be  checked 
in  the  more  meaningful  environment  of  its  intended  appli- 
cation. An  accuracy  check  is  built  right  into  the  component 
matrix  approach  to  system  identification:  all  of  the  com- 
ponent matrices  of  the  pseudo  modes  should  be  null.  Also, 
the  computational  accuracy  should  be  tested  on  computers 
with  shorter  word  length  than  possessed  by  the  CDC  6600 
computer  used  for  these  tests.  Computers  used  for  on-line 
control  applications  probably  wouldn't  have  word  lengths 
with  48  bits  of  data  storage  like  the  CDC  6600  has. 
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The  following  procedures  are  therefore  recommended  for 


VANINV  itself 


1.  Test  VANINV  with  repeated  and/or  complex  eigen- 
values to  determine  their  effect  on  computation 
time  and  accuracy. 

2 . Determine  where  most  of  the  computation  time  is 
used  on  VANINV  and  improve  the  efficiency  where 
practical. 

3.  Eliminate  the  calls  to  subroutines  that  are  made 
only  one  or  two  places  in  VANINV  and  replace  them 
with  the  coded  subroutines,  since  all  the  sub- 
routines have  been  debugged. 

4.  Streamline  VANINV  by  using  vector  array  addressing 
to  reduce  array  index  computation  time. 

5.  Incorporate  a row  jumping  capability  into  VANINV 
to  eliminate  duplicated  calculations  for  repeated 
and/or  complex  eigenvalues.  This  is  illustrated 
in  Figure  8,  which  shows  a very  common  occurrence 
for  the  component  matrix  approach.  Every  complex 
conjugate  pair  of  eigenvalues  in  the  original 
system  will  be  repeated  in  the  2n  x 2n  Vandermonde 
matrix  of  the  augmented  system.  The  resulting 
calculations  that  VANINV  would  perform  for  these 
eigenvalues  is  shown  by  Figure  8.  Note  that  the 
lower  half  of  the  matrix  shown  in  Figure  8 is 
simply  the  complex  conjugate  of  the  upper  half. 
Thus  considerable  computation  time  could  be  saved 


mm 


t 

Hi 


by  incorporating  a row  jumper  which  would  simply 
store  the  complex  conjugate  of 

Cil  _ d N£CS) 

WU  ' 3i  UTTsT 

xl 

directly  into 

C2)  _ d N£(-S) 

Wl£  ds  D0(S) 

*2 

and  thus  eliminate  half  the  calculations.  Also 
note  that  a repeated  eigenvalue  in  the  original 
system  would  occur  four  times  in  the  2n  x 2n  Vander- 
monde matrix,  as  Figure  9 shows.  Each  column  is 
composed  of  identical  elements  except  for  the 
number  of  derivatives  that  must  be  taken  for 
evaluation.  If  a row  jumper  with  some  memory  were 
added  to  the  routine,  successive  derivatives  could 
be  taken  to  evaluate  the  elements  of  each  column. 
This  would  require  only  1+1+1  = 3 derivatives  for 
each  column  shown  rather  than  3+2+1  * 6 derivatives 
as  VANINV  must  do  presently.  (Without  the  row 
jumper,  VANINV  simply  evaluates  each  element  as 
it  comes  to  it,  without  regard  for  any  calculations 
done  previously.) 

In  addition  to  the  recommended  procedures  for  VANINV 
described  above,  two  comments  are  made  regarding  the  com- 
ponent matrix  approach  to  system  identification: 
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1.  Many  of  the  matrices  resulting  from  the  com- 
ponent matrix  approach  are  very  sparse  (eg:  Eqs 
(D2-13)  , (D2-14) ) , and  thus  the  method  lends 
itself  readily  to  sparse  matrix  operations  such 
as  those  available  with  SOFE  (Ref  7) . 

2.  The  calculation  structure  of  both  the  component 
matrix  approach  and  VANINV  lend  themselves  to 
parallel  calculations.  An  implementation  of 
these  methods  on  a parallel  architecture  com- 
puter has  the  potential  for  being  very  fast  in 
real  time. 

Thus  it  may  be  seen  that  a refinement  of  VANINV  may 
well  have  a place  in  future  system  identification  processes. 
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APPENDIX  A 


COMPUTER  PROGRAMS  FOR  CSAKI'S 
GENERALIZED  INVERSE  VANDERMONDE 


ALGORITHM 


Preface  for  Appendix  A 


This  appendix  contains  the  subroutines  written  to  implement  the 
generalized  method  for  finding  an  inverse  Vandermonde  matrix  that  is 


serve  as  a combination  user/programmer  guide  for  further  investigation 


of  the  proposed  alternate  method  of  system  identification  which  is 


outlined  in  the  main  text 


Subroutines  Called 


t Not  used  by  any  of  the  subroutines  listed,  but  pertinent  to  the 
use  of  VANXNV. 

t Function  subprogram. 


h eu 


A-5 


A- 7 


A-10 
A- 11 


A-13 


A-15 


A-18 


A-19 


A-21 


A- 2 3 
A -2  4 


A- 26 


A-28 


A-30 


A-32 


A- 34 


A- 36 


A-40 
thru 
A -4  6 


A- 3 


Subroutine  CGAIN 


Subroutine  CGAIN  multiplies  each  coefficient  of  a polynomial  by 
specified  gain.  The  coefficients  and  the  gain  may  be  complex. 

Subroutine  statement:  SUBROUTINE  CGAIN  (AR.AI,NSA,AKR,AKI) 
Subroutines  called:  None 

Arrays  containing  the  real  and  imaginary  parts  of 
the  polynomial  coefficients 

Number  of  storage  locations  required  for  the 
polynomial  coefficients 

Real  and  imaginary  parts  of  the  polynomial  gain 


Variables: 

ARl 

All 

NSA  - 

AKRl 

AKI1 


A-4 


1 ’ • +"• 


- V V , s 


■■ 


■■■■b l 


^Subroutine  CGAIN  (AR,AI ,NSA,AKR, AKI)  J 

Is  the  gain  real? 


' Close 

>|aki| 

v ? . 


No ’ . . . Multiply  by 
complex  gain. 


Is  the  GAIN 
unity? 


AKR  = 1 
\ ? / 


TR  * AR( I) 

TI  - AI (I) 

AR(I) 
AI  (I) 

■ TR*AKR-TI*AKI 
- TR*AKI+TI*AKR 

No  ...  Multiply 
by  real  GAIN. 


Done? 


I - 1+1 


I 1NSA  ^ 


I)  * AR( I) *AKR 
I)  - AI (I) *AKR 


Done? 


I > NSA 

V' 


26  ^ Update  GAIN. 

AKR  - 1. 

AKI  - 0. 


Figure  A-l.  Flowchart  for  subroutine  CGAIN. 
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Subroutine  COPYPOL 


Subroutine  COPYPOL  copies  the  coefficients  of  the  input  polynomial 


into  the  storage  locations  for  the  coefficients  of  the  output  polynomial 


Also,  the  input  gain  is  copied  into  the  output  gain  storage.  The  coef- 


ficients and  gains  may  be  complex. 


Subroutine  statements  SUBROUTINE  COPYPOL  CPIR,PII,GIR,GII,N,POR,POI, 

GOR, GOII 


Subroutines  called:  None 


Variables: 


PIR) 

PHI  Arrays  containing  the  real  and  imaginary  parts  of  the 

PORI  * coefficients  of  the  input  (PLxl  and  output  CPOx)  polynomials 

POI1 


N - Number  of  storage  locations  needed  for  the  polynomial 

coefficients.  (Note  that  N is  one  integer  value  larger 
than  the  polynomial  order.) 


Reed,  and  imaginary  parts  of  the  gedns  of  the  input  CGIx] 
and  output  (GQxl  polynomials 


P • 


■ 


* 


Subroutine  ESORT 

Subroutine  ESORT  sorts  a collection  of  real  and  complex  eigenvalues 
according  to  the  following  specifications: 

1.  Put  in  order  of  decreasing  complexness;  i.e.,  complex  before 
real. 

2.  Order  each  section  in  decreasing  order  of  multiplicity. 

3.  Order  complex  conjugates  with  the  eigenvalue  having  the 
positive  imaginary  part. 

4.  For  equal  multiplicities,  order  in  terms  of  decreasing 
magnitude  of  the  real  part;  or  if  equal,  place  eigenvalues 
with  positive  real  parts  before  those  with  negative  real 
parts,  and  order  in  term®  of  the  imaginary  part. 

Subroutine  ESQRT  sweeps  through  the  collection  and  compares  each  pair 

of  adjacent  eigenvalues  to  determine  whether  or  not  they  need  to  be 

swapped  to  get  them  in  the  proper  relative  order.  Each  time  a swap  is 

performed,  a "swap  counter"  is  incremented.  At  the  end  of  each  sweep 

through  the  eigenvalue  array,  the  swap  counter  is  tested  to  see  if  any 


swaps  were  made  on  that  sweep.  If  there  were  swaps,  the  counter  is 
reset  to  zero  and  another  sweep  through  the  array  is  started.  Zero 
swaps  during  any  sweep  through  the  array  indicates  that  all  the  eigen- 
values are  in  the  proper  order,  so  the  subroutine  returns  the  rearranged 
eigenvalue  array. 


Subroutine  statement:  SUBROUTINE  ESORT  CEIGR.EIGI ,KI,M) 

Subroutines  called:  None 

Variables: 


EIGR)  ^ Arrays  containing  the  real  and  imaginary  parts  of 
EIGI)  the  collection  of  eigenvalues 

KI  ■ Array  containing  the  corresponding  eigenvalue 
multiplicities 


A- 8 


I 


M * Number  of  different  eigenvalues  in  the  array  (complex 
conjugates  are  considered  to  be  different) 


MDIS  ■ "Message  Disable"  flag,  passed  in  blank  common.  ESORT 

normally  prints  an  informative  diagnostic  if  two  identical 
eigenvalues  are  encountered.  If  the  user  doesn't  want  the 
message  printed,  setting  MDIS  equal  to  any  positive  integer 
will  disable  the  printer.  Setting  MDIS  to  zero  or  a 
negative  integer  enables  the  printer. 


<• 


routine  ESORT  (EIGR,EIGI,KI ,Ml 


^ I Zero  the  Swap  Counter 


I SWAP =0 


Point  to  the  I'th  Eigenvalue 


Point  to  the  next  Eigenvalue 


IN-I+1 


Vre  there  Two  Eigenvalues 
to  Compare? 


Test  I'th  and  In'th 
Eigenvalues  to 
determine  if  they 
should  be  swapped. ^ 


Any  Swaps  performed 
10  w this  Sweep? 

ISWAP  >\vYes . • .Make  another  Sweep 

>0  ? 

V ? / 


Return 


EIGRT-EIGRCI) 

EIGIT-EIGI(I) 

KIT-KI(I) 

EIGR(I) =EIGR(IN) 
EIGI(I)-EIGKIN) 
KI  (I)-KI (IN) 
EIGR(IN) “EIGRT 
EIGI (IN) “EIGIT 
KI(IN)-KIT 

♦ ' 

|ISWAP-ISWAP+1 


140  l Look  at  Next  Eigenvalue  Pair 


I-I+l 


*See  page  A-ll  for  flowchart. 


Figure  A-3 . Flowchart  for  subroutine  ESORT. 


Bogin  ESORT 
Toot  Logic 


Subroutine  EVAL 


Subroutine  EVAL  finds  the  value  of  the  input  polynomial  at  the 
specified  value  of  the  variable.  The  input  polynomial  coefficients  and 
the  value  of  the  variable  may  be  complex. 


Subroutine  statements  SUBROUTINE  EVAL  (AR,AI,AKR, AKI,EIGR,EIGI, 

NS ,BR,BI) 


Subroutines  called: 


Variables: 


CGAIN 


AR)  m Arrays  containing  the  real  and  imaginary  parts 
AI)  of  the  coefficients  of  the  input  polynomial 

AKR)  m Real  and  imaginary  parts  of  the  gain  of  the  input 
AKI)  polynomial 


EIGR)  _ Real  and  imaginary  parts  of  the  value  at  which  the  poly- 
EIGI)  nomial  is  to  be  evaluated  (must  be  chosen  by  the 
calling  routine) 

NS  ■ Number  of  storage  locations  required  for  the  input 
polynomial;  NS*N+1,  where  N is  the  order  of  the 
polynomial 

BR) 

“ Real  and  imaginary  parts  of  the  returned  value 


Find  Polynomial  order. 


N - NS- l 


Put  conatant  term  in  output  storage. 


BI-AI(NS) 


Check  for  more  terms . 


Nora  terms  found  . . . Continue 
evaluation. 


Costplex  eigenvalue? 


Yes  ...  Use  complex 
evaluation . 


No  ...  Use  real 
evaluation- 


| EI-EIGI 


Add  next  higher 
order  term  to  output. 


MUi-AJM.N)| 

UN-AI(N) 


Add  next  higher 
order  term  to  output. 


BR-BR*AR(N) *ER 
BI-BI*AI(N) »ER 


BI-BI*AJWEI*AIN*ER 


ate  evaluation  terms. 


Update  evaluation 
term. 


rR-ER*EIGR-EI«lIC 

ri><R*cxGi«c:*txa 

ER-TR 

II**! 


I £R-tR*EISR 


Any  more  terms7 


Any  more  terms? 


J • J*l 


Multiply  by  complex  gain 


Figure  A-5.  Flowchart  for  subroutine  EVAL. 
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Subroutine  EXPAN DC 


Subroutine  EXPANDC  finds  the  coefficients  for  the  polynomial  expan- 
sion of  a collection  of  root  factors.  EXPANDC  reads  the  roots  from  a 
root  array  and  returns  the  coefficients  in  a polynomial  coefficient 
array.  The  coefficient  for  the  highest  power  of  the  variable  is  the  first 
element  in  the  coefficient  array.  The  routine  handles  both  real  and 
complex  roots  and  does  not  require  that  comp lex  roots  have  conjugate 
pairs.  The  maximum  size  of  the  polynomial  found  is  limited  only  by  the 
storage  allocated  for  CR  and  Cl  in  the  calling  program. 


Subroutine  statement:  SUBROUTINE  EXPAN DC  (ROOTR, ROOTI ,NS,BRf BI,CR,CI, 

CKR, CKI) 


Subroutines  called:  POLYMC 


Variables: 

ROOTR) 

ROOTI) 


NS 


BR) 

BT) 


Arrays  containing  the  real  and  imaginary  parts  of  the 
roots  whose  factors  are  to  be  expanded  into  a 
polynomial 

The  number  of  storage  locations  required  for  the  output 
polynomial;  equal  to  NF+1  where  NP  is  the  number  of 
factors 

NS-dimensioned  temporary  storage  for  intermediate 
results  of  the  polynomial  expansion 


CR)  k NS-dimensioned  arrays  containing  the  real  and  imaginary 
CZ)  parts  of  the  coefficients  of  the  output  polynomial 


CKR)  m Real  and  imaginary  parts  of  the  output  polynomial  gain 
CKI)  multiplier;  always  equal  to  1.  and  0.,  respectively 

CLOSE  “ User  specified  value  used  as  a null  tester.  If  the 

magnitude  of  any  value  tested  is  less  than  CLOSE,  that 
value  is  set  to  zero.  CLOSE  is  stored  in  blank 
common  and  must  be  specified  in  the  main  program. 
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i LMI-n  ifii 

nrodro. 


:«*t  Sac  restart  ■ 

It  UMf«  «c*  no  f ACT  at » . ;tint  trrat  miiw«  mA  itsp. 
tf  :n*r*  ;c  an..-  an*  factor.  return  ^B*4i.«c*ly. 

1 :S  awrt  ac*  cko  ar  aor*  factor*,  cwuau*  <iU  *.<pansion. 


»nnt  trtot 


*a»* 

xjc:«i 


Tt  JMuury  pact  of  um  coot  iifBUicMtf 
tot  up  MOUr  accordingly. 


tot  ...  tot  up 

MOLT.  ignoring  tto  toll  racy 


Figure  A-6.  Flowchart  for  subroutine  EXPANDC 


Subroutine  FRACDIF 


Subroutine  FRACDIF  accepts  a polynomial  fraction  and  finds  its 
derivative  using  the  quotient  rule: 

dv  du 
d ,u.  u dx  V lx 


The  coefficients  of  the  input  numerator  and/or  denominator  polynomial 
may  be  complex. 


Subroutine  statement:  FRACDIF  (PUR,PUI,PVR,PVI ,NUS,NVS,PNR,PNI,PDR, 

PDI,NNS,NDS,NT,UR,UI,VR, VI,RNR,RNI,DR,DI, 
TEMP1R, TEMP1I ,TEMP2R, TEMP2I , TEMP 3R, TEMP3 I) 

Subroutines  called:  POLYDIF,  POLYMC,  POLY SUB 


Variables: 

PUR,PUI)  Arrays  containing  the  real  and  imaginary  parts  of  the 
PVR* P VI)  coefficients  of  the  input  polynomials,  U and  V 


NUS) 

NVS) 

PNR/PNI) 

PDR,PDI) 


NNS) 

NDS) 


NT 


Number  of  storage  locations  needed  to  store  the 
coefficients  of  the  input  polynomials 

Arrays  containing  the  real  and  imaginary  parts  of 
the  coefficients  of  the  output  numerator  polynomial 
(PNx)  and  denominator  polynomial  (PDx) 

Number  of  storage  locations  needed  to  store  the 
coefficients  of  the  output  numerator  find  denominator 
polynomials,  where: 

NNS  * NUS  + NVS 

NDS  - 2NVS  - 1 

Number  of  temporary  storage  locations  needed  for  the 
numerator  calculations 

NT  - NON  (L)  + NOD  (I) 

where:  NON (L)  is  storage  required  for  L'th  numerator 

polynomial 

NOD (I)  is  storage  required  for  fth  demoninator 
polynomial 
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Real  and  .imaginary  parts  of  the  gains  associated  with 
polynomials  U,  V,  PNx,  PDx 


UR#  UI) 

VR,  VI) 

FUR,  RNI) 

OR# 01} 

TEMP1R,TEMP1I1 

TEMP 2R, TEMP 21)  - Temporary  storage  locations  for  numerator  calculations 

TEMP3R,TEMP3I) 

, 


Subroutine  FRACDIF IPUR, PUI , PVR, PVI ,NUS , NVS ,PNR, PNI ,PDR,PDI , 
NNS,NDS,NT,UR,UI,VR,VI,RNR,RNI,DR,DI, 
TEMP1R, TEMPI I , TEMP2R, TEMP 2 I , TEMP 3R, TEMP 3 I). 


Figure  A-7.  Flowchart  for  subroutine  FRACDIF. 
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Subroutine  QUT2C 

Subroutine  0UT2C  prints  out  the  elements  of  a two-dimensional 
array  and  labels  each  element  with  the  array  name  and  its  position  in 
the  array.  The  real  and  imaginary  parts  of  the  complex  array  are 
treated  separately. 

Subroutine  statement:  SUBROUTINE  0UT2C  (NR ,NC ,AR, AI , IDR, IDI) 
Subroutines  called:  None  , 

Variables: 

NR  * Number  of  rows  to  be  printed  out 

NC  = Number  of  columns  to  be  printed  out 

AR)  m Arrays  containing  the  real  and  imaginary  parts  of 
All  the  array  to  be  printed 

IDR}  Hollerith  labels  for  the  real  and  imaginary  parts  of  each 
IDI]  * element  printed;  these  must  be  specified  in  Hollerith 
format  by  the  calling  program 


20 


Subroutine  POLYADD 


— mu— m 


Subroutine  POLYADD  adds  two  polynomials,  POLYA  and  POLYB , and 
stores  the  result  in  POLYC.  The  coefficients  of  the  polynomials  may 
be  complex.  After  the  addition  is  complete,  the  first  coefficient  of 
the  output  is  tested  to  determine  if  it  is  unity.  If  it  is  not,  subroutine 
UNITY  is  called  to  unitize  it. 


Subroutine  statements  SUBROUTINE  POLYADDCAR,AI,BR,BI,CR,CI,NSA,NSB,NSC, 

AKR, AKI , BKR, BKI , CKR, CKI ) 


Subroutines  called: 


CGAIN,  UNITY 


Variables: 


AR1 

All 

BRJ  Arrays  containing  the  real  and  imaginary  parts  of 

BI)  " the  coefficients  of  POLYA,  POLYB,  and  POLYC 

CR1 

oil 


NSAl 

NSB1 

NSC). 


Number  of  storage  locations  required  for  the  coefficients 
of  POLYA,  POLYB,  and  POLYC 


AKRL 

AKI) 

BKR1 

BKI) 

CKR) 

CKIl 


Real  and  imaginary  parts  of  the  gain  terms  for 
POLYA,  POLYB,  and  POLYC 
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Subroutine  POLYADD  (AR, AI , BR , BI , CR, Cl , NSA , NSB , NSC , 
AKR, AKI ,BKR,BKI ,CKR,  CKI).  . 


y Set  output  gain. 

CKR*  1. 

CKI » Q. 

^ Multiply  POLYA  by  its  gain  if  not  Unity. 


Call  CGAIN (AR,AI , NSA, AKR, AKI) 


Multiply  POLYB  by  its  gain  if  not  Unity. 


Call  CGAIN (.BR,BI,NSB,BKR,BKI) 


d#  Zero  the  output  array. 


CR(_I)  =Q. 
CIU1-Q. 


1=  1+1 


Done? 


I > NSC 


Add  POLYA  to  the 
Output  Array 

y ::  zz 

Add  POLYB  to  the 
Output  Array 


I If  first  coefficient  of  output  is  not 
I Unity,  unitize  it. 

JL Real  part  equal  to  one? 

CRCllV  N° 

-1.  > 1 


i Yes  . . . Imaginary  part  equa 
yS\to  zero? 

CICllV  No 

-o.  s a 


[Call  UNITY (CR,CI ,NSC,CKR , CKI) 


Return 


See  page  A-24  for  flowcharts . 
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Add  Polynomials  to  I 
the  Output  Array  I 
POLYADDJ^J 

Add  POLYA  to  the  Output  Array. 

Map  first  coefficient  of  POLYA  to  proper 
^ f location  in  Output  Array. 

NSTART  »N  5C-N  SA- 1 ~| 

1 Point  at  the  first  coefficient  of  POLYA. 


Add  each  coefficient  of  POLYA  to  its  proper 
Location  in  the  Output  Array. 


CR(J)»AR(I) 
CI(J)  «AI(I] 


Done  adding? 


Yes  ...  Add  POLYB  to  the  Output  Array. 

Map  first  coefficient  of  POLYB  to  proper 
location  in  Output  Array. 


NSTART  « NSC-NSB+1 


Point  at  the  first  coefficient  of  POLYB. 


Add  each  coefficient  of  POLYB  to  its  proper 
location  in  the  Output  Array. 


CI(J)-CI(J)+BI(I) 


Done  adding? 


J > NSC 


End  of  Adding  Poly- 
nomials  to  the 
Output  Array 


Figure  A-ll.  Flowchart  for  addition  steps  of  subroutine  POLYADD . 
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Subroutine  POLYDIF 


Subroutine  POLYDIF  accepts  a polynomial  input  and  returns  the 
derivative.  The  coefficients  of  the  input  polynomial  may  be  complex. 


Subroutine  statement:  SUBROUTINE  POLYDIF (POLYR.POLYI, NS, DPOLYR, 

DPOLYI , N , PKR , PKI 1 

Subroutines  called:  CGAIN 

Variables: 

POLYR) 

POLY I) 

DPOLYR1 
DPOLYI) 

NS) 

N ) 


PKR) 

PKI) 


Arrays  containing  the  real  and  imaginary  parts  of  the 
coefficients  of  the  input  polynomial 

Arrays  containing  the  real  and  imaginary  parts  of  the 
coefficients  of  the  output  (derivati’- - polynomial 

Number  of  storage  locations  required  for  the 
coefficients  of  POLY  and  DPOLY  (Note  that: 

NS  =*  NI+1,  N = NI}  where  NI  is  the  order  of  the 
input  polynomial). 

Real  and  imaginary  parts  of  the  polynomial  gain 
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Subroutine  POLYMC 

Subroutine  POLYMC  finds  the  coefficients  of  the  product  of  two 
polynomials.  The  routine  handles  all  combinations  of  real  and  complex 
input  coefficients.  The  order  of  the  polynomials  involved  is  limited 
only  by  the  storage  allocated  for  them  by  the  calling  program. 


Subroutine  statement: 


SUBROUTINE  POLYMC (AR,AI,BR,BI,CR, Cl, NSA,NSB, NSC, 
AKR , AKI , BKR , BK I , CKB , CKI 1 


Subroutines  called: 
Variables : 


UNITY 


AR) 

All 


Arrays  containing  the  real  and  imaginary  parts  of  the  coef- 
ficients of  the  first  polynomial  multiplicand,  POLYA 


BR1  m Arrays  containing  the  real  and  imaginary  parts  of  the  coef- 
BI]  * ficients  of  the  second  polynomial  multiplicand,  POLYB 


CR] 

CI1 


Arrays  containing  the  real  and  imaginary  parts  of  the  coef- 
ficients of  the  polynomial  product,  POLYC 


NSA1 

NSB] 

NSC1 


Number  of  storage  locations  required  for  the  coefficients 
of  POLYA,  POLYB,  POLYC;  each  is  one  integer  larger  than 
the  corresponding  polynomial  order  (Note:  NSC  ” NSA+NSB-1) 


AKR) 

AKI1 

BKR] 

BKI) 

CKR) 

CXI1 


Real  and  imaginary  parts  of  the  gains  associated 
with  POLYA,  POLYB,  POLYC 
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M,  ^ * 


of  the  gain  of  POLYB  and  calling  subroutine  POLYADD.  The  coefficients 


SUBROUTINE  POLYSUB (AR,AI ,BR,BI ,CR,CI,NSA,NSB,NSC 
AKR,AKI,BKR,BKI, CKR,CKI} 


POLYADD 


Arrays  containing  the  real  and  imaginary  parts  of  the 
coefficients  of  the  polynomials  operated  upon  according 
to  POLYOPOLYA-POLYB 


Numher  of  storage  locations  required  for  the  coefficients 
of  POLYA,  POLYB,  and  POLYC,  respectively.  Note  that  NSC 
equals  the  larger  of  the  two  values  NSA  and  NSB. 


AKR) 

AKI} 

BKRj  Real  and  imaginary  parts  of  the  gain  multipliers  for 
BKIl  * POLYA,  POLYB,  and  POLYC,  respectively. 

CKR) 

CXI) 


Subroutine  RAT 


Subroutine  RAT  rationalizes  a complex  fraction  by  multiplying  and 
dividing  the  fraction  by  the  complex  conjugate  of  the  denominator . The 
denominator  coefficients  are  tested  for  nullity  prior  to  the  last  step 
of  rationalization,  and  if  they  are  null,  an  error  message  is  printed 
and  execution  is  terminated. 


Subroutine  statement:  SUBROUTINE  RATCRNR,RNX,DR,DI,RR,RI) 
Subroutines  called:  None 


Variables : 


RNR) 

RNI) 


Real  and  imaginary  parts  of  the  input  numerator 


DR} 

DI} 


Real  and  imaginary  parts  of  the  input  denominator 


RR} 

RI1 


Real  and  imaginary  parts  of  the  rationalized  output 
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Subroutine  ROOTAY 


Subroutine  ROOTAY  expands  a collection  of  eigenvalues  with  their 
associated  multiplicities  into  an  array  containing  one  root  for  each 
occurrence  of  each  eigenvalue.  Both  real  and  complex  eigenvalues  may 
be  used.  The  routine  forms  the  root  array  beginning  with  the  IB'th 
eigenvalue  in  the  input  array  and  uses  all  the  remaining  eigenvalues 
in  the  array.  Thus  the  input  array  must  have  the  eigenvalues  arranged 
so  that  any  eigenvalues  to  be  excluded  from  the  expansion  precede  the 
IB'th  eigenvalue. 


Subroutine  statement:  SUBROUTINE  ROOTAY (EIGR,EIGI,K,M, IS, IB, 

ROOTR, ROOTI) 


Subroutines  called: 
Variables: 


None 


EZGR1 

EIGI) 


Arrays  containing  the  real  and  imaginary  parts  of  the 
input  collection  of  eigenvalues 


K ■ An  array  containing  the  multiplicities  of  the  eigenvalues 
M * The  number  of  different  eigenvalues 
IS  - 


The  number  of  storage  locations  needed  to  contain  the 
resulting  array  of  roots; 


note  that  IS 


£ • 


the  summation  of  the  eigenvalue 


i-1 


multiplicities 


ROOTR)  Arrays  containing  the  real  and  imaginary  parts  of  the 
ROOTI)  roots  in  the  output  array 


IB 


The  number  of  the  eigenvalue  with  which  to  begin  forming 
the  root  array. 
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Subroutine  ROOTAY (EIGR,EIGI ,K,M,IS, IB,ROOTR,ROOTI) 


J - J+l 


IE-IE+1 


Point  at  the  IR'th  location  of  the 
root  array. 


Start  with  the  IB'th  eigenvalue. 


Determine  the  multiplicity  of  the 
IE ' th  eigenvalue . 


Copy  the  eigenvalue  into  the 
root  array  L times. 


R00TR(.IR1  -EIGRUE1 
ROOTI  (IRi-EIGI  CIE1 


IR-IR+1 


Done  with  root? 


J > L 


Done  with  array? 


IE  _>  M 


Return 


Figure  A- 16.  Flowchart  for  subroutine  ROOTAY. 
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Subroutine  UNITY 


. 


Subroutine  UNITY  operates  on  the  input  polynomial  to  unitize  the 
coefficient  of  the  highest  order  variable.  This  is  accomplished  by 
extracting  the  coefficient  of  the  highest  order  term  as  a gain  multi- 
plier. The  polynomial  gain  is  then  updated  to  reflect  the  extraction. 
The  input  coefficients  may  he  real  or  complex.  In  the  event  that  the 
first  coefficient  of  the  input  polynomial  is  zero,  it  is  ignored  and 
the  next  one  is  tested.  If  all  the  coefficients  are  zero,  the  gain 
and  constant  term  are  set  to  zero,  and  a warning  message  is  printed 
if  the  printer  has  not  been  disabled. 


Subroutine  statement:  SUBROUTINE  UNITY (AR, AI, NS, AKR,AKIl 
Subroutines  called:  None 


Variables: 


AR]  m Arrays  containing  the  real  and  imaginary  parts  of  the 
AI]  coefficients  of  the  input  polynomial 

NS  * Number  of  storage  locations  required  for  the  coefficients 
of  the  input  polynomial 


AKR] 

AKI) 


Real  and  imaginary  parts  of  the  polynomial  gain 


CLOSE  * User  specified  magnitude  used  in  testing  the  significance 
of  variables.  If  the  magnitude  of  the  tested  variable  is 
less  than  CLOSE,  it  is  assumed  to  be  zero.  CLOSE  is  stored 
in  blank  cannon. 

MDIS  ■ User  specified  printer  disable  flag.  If  MDIS  is  greater 
than  zero  all  print  steps  using  the  flag  are  disabled. 

MDIS  is  stored  in  blank  common. 
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Subroutine  VAN IN V 

Subroutine  VANINV  is  an  implementation  of  Csaki's  generalized 
algorithm  for  finding  the  inverse  of  the  Vandermonde  matrix.  It  accepts 
an  ordered  array  of  real  and  complex  eigenvalues  and  their  multiplicities, 
and  returns  the  inverse  Vandermonde  matrix.  Prior  to  calling  VANINV,  the 
eigenvalues  must  be  arranged  in  the  order  described  for  subroutine  ESORT. 


Subroutine  statement: 


SUBROUTINE  VANINV CEIGR,EIGI,KI,N,M, NS, MS, WR,WI, 
NW1,NW2 ,AR,AI,FR,FI,NR,NI,DR,DI,GR,GI,NOD,NON, 


PNR,PNI ,PDR,PDI,PUR,PUI ,PVR,PVI ,TlR,TlI ,T2R,T2I , 
T3R,T3I,DKR,DKI1 

Subroutines  called:  COPYPOL , EVAL , EXPANDC , FRACDIF , RAT , ROOTAY 

Function  subprograms  called;  NFACT 


Variables: 


11  Formed  parameters 


EIGR) 

EIGI) 

KI 

N 

M 


Arrays  containing  the  real  and  imaginary  parts  of  the  M 
different  eigenvalues  of  the  N'th  order  system 


Array  containing  the  multiplicities  corresponding  to  the 
system  eigenvalues 


System  order  {note  that  N 


KU1 


i*l 


Number  of  different  eigenvalues  (complex  conjugates  are 
considered  to  be  different) 


NS)  ^ Storage  size  parameters,  equal  to  N+l  and  M+l, 

MS)  respectively 

WR)  N x n arrays  into  which  the  real  and  imaginary  parts  of 
WI)  = the  inverse  Vandermonde  matrix  are  placed  by  VANINV 

NWl]  = Storage  size  parameters,  equal  to 

NW2)  * (2** (KI  -1))*(N-KI  )+l  and  2N,  respectively 

max  max 

AR)  N-dimensioned  arrays  used  for  intermediate  polynomial 

AI)  * calculations  (short  form  of  P0LYA1 
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FR) 

FI)  NS-dimensioned  arrays  used  for  intermediate  polynomial 
GR)  calculations  (short  forms  for  POLYF  and  POLYG) 

Gil 

NR)  _ (N  x N) -dimensioned  arrays  used  to  store  the  real  and 
NI)  imaginary  parts  of  the  N truncated  numerator  polynomials 

DR)  _ (M  x MS) -dimensioned  arrays  used  to  store  the  real  and 
DI)  imaginary  parts  of  the  M denominator  polynomials 

NOD  * M-dimensioned  array  used  to  store  the  polynomial  order 
associated  with  each  of  the  M denominator  polynomials 

NON  = N-dimensioned  array  used  to  store  the  polynomial  order 
associated  with  each  of  the  N numerator  polynomials 


PNR,  PNI) 
PDR, PDIl 
PUR, PUI1 
PVR,PVI) 


T1R,T1TJ 

T2R,T2I) 

T3R.T3I1 


AKR.AKI) 
FKR, FKI) 
GXR,GKI) 
UR,UI  1 
VR, VI  1 
RNR, RNI) 
RDR.RDI) 


GNLRl 

GNLI) 


NDIV 


NDIVS 

RC 


working  storage  for  polynomials  used  in  differentiating 
a polynomial  fraction  by  means  of  the  quotient  rule: 


dv 


— (-) 
dx  v 


U dx  “ V ^dx1  POLYN 


POLYD 


Note  that  all  of  these  must  be  dimensioned  at  least  NWl 


NW2-dimensioned  arrays  used  for  temporary  storage 


Real  and  imaginary  parts  of  the  gains  associated  with 
POLYA,  POLYF,  POLYG,  POLYU,  POLYV,  POLYN,  POLYD 


Real  and  imaginary  parts  of  the  gains  associated  with 
each  of  the  M denominator  polynomials.  Since  in  all 
cases:  GNLR?*1 . and  <3JLI=0.,  they  are  defined  with  a 
data  statement. 


Number  of  derivatives  that  must  be  taken  of  the  current 
inverse  Vandermonde  element,  W(J,L) , before  evaluation. 


Storage  for  NDIV  so  initial  value  can  be  recalled 


Row  coefficient  for  J'th  row  of  I'th  block  of  the  inverse 
Vandermonde  matrix 

Note  that  RC  - (K1CI)  jj  i (NDIVS)  '. 


A- 38 


v *»<•*■* 


NVS1 
NNS)  - 
NDS) 
NT) 


Dummy  varieties  used  to  pass  between  subroutines  the  cur- 
rent size  of  the  storage  locations  of  POLYN  and  POLYD 


RNUMR, RNUMI)  = Real  and  imaginary  parts  of  the  polynomial  fraction  for  an 
DENR,  DENI  ) element  of  the  inverse  Vandermonde  prior  to  rationalization 


ZR)  _ Real  and  imaginary  parts  of  an  element  of  the  inverse 
ZI)  Vandermonde  matrix  after  rationalization  of  the  fraction 
but  prior  to  multiplication  by  the  row  coefficient 
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Subroutine  VANINV(EIGR,EIGI,KI,N,M,NS,MS,WR,WI,NW1,NW2,AR,AI,FR,FI, > 
NR, N I , DR, DI , GR, GI ,NOD , NON , PNR , PNI , PDR.PDI , PUR, PUI , 
PVR,PVI,T1R,T1I,T2R,T2I,T3R,T3I,DKR,DKI) 


Check  for  trivial  solution  (only 
one  eigenvalue) . 


"m  - 1 
< 0 


Trivial. 


No  ...  Non 
trivial . 


iSet  output 
values . 
WR(1,1)«1.  I 
WI(1,11*0.  I 


Sac  up  N Denominator 
Polynomials  for  Subroutine 
VMtNVi  DSC)  .DIUl.NOOUl 
_ (P/0  VANINV) 


Set  up  root  errey  Cot  I'th  polynomial . 


11-2 

IT-W-Kim 


11  POOT1Y (SIGP.EIGI  .HI . H,  IT.  II.  11, All  1) 

V finJ  PCI)  polynomial  froe  root  errey. 


ceil  txtMiPc<M>.Ai.if.ni.ri.<a>.Gi.gq>.gci)  II 

W store  returned  poiynonisl  ee  D(I) i end  ite  9eln. 

I**  1 1 


lOlU.D-OllR) 


Tee  ...  Store 


|DKiii)e<aa  | 

t Store  the  order  of  the  I'th  polynomial. 


ate  eigenvalues  to  find  neat  polynoeiel. 


IISIT-CISIU) 
KT  • KX(1) 


tlOI(K)  -tICKI 
KI  (K)  -K1  ( L) 


Done  rotating? 


cxeiaixisrr 

kkd-ci 


Done  finding  denominator  polynomials? 


End  setup  of 
Denominator 
Polynomials 


Figure  A-19.  Flowchart  for  denominator  polynomial 
setup  steps  for  subroutine  VANINV. 
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Setup  N Numerator 
Polynomials  for  Subroutine 
VANINV;  NR(I) ,NI(I) ,NON(I) 

(P/0  VANINV)  _ 

1 Find  characteristics  equation , 


Call  ROOTAY (EIGR,EIGI ,KI ,M,N, 1 , AR, AI) 


Call  EXPANDC ( AR , AI , NS , FR , FI , GR , GI , GKR , GKI ) 


w Begin  polynomial  truncation. 


[nl 

.=N 

* 

I L» 

i] 

j 

f 

j 

■n 

J » J+l 


NR(L , J) =GR(J) 

NI (L, J) *0 . 

g7  y Done  truncating  L'th  polynomial? 

0 /\ 

No>^  J _>  NL  \ 

T Yes  ...  Store  polynomial  order. 
NON (L) *NL 

* 

NL  - NL-1 

a W Done  with  truncations? 

— /l  > N \ 


Figure  A- 20. 


End  Setup  of 
Numerator 
Polynomials^ 


Flowchart  of  numerator  polynomial  setup 
steps  for  subroutine  VANINV. 


■pacify  Starapa  Required  for  Evaluation  ao 
rropar  tolyrwai  il  a for  maarator  and 
ulna  tor  into  wort  lop  krraya.  * 


| tl-Eiaill)  | 

^ rraluata 

I call  rviu.tnnt.»m.i>ot.mi,i«.ii.»i«, 

i Cvaluata 

Tag  tvmcm«.ri>i.tp«.tci.t»,fx.»o«. 


row  eoafficiant  to  pat  'V  tar 


ti(n,u»ani 


Taa  ...  Daaa  with 


with  natxix? 


t |M  papa  MS  (or  flowchart. 


td  Solvinp  for  Zlawantai  W»(J.L)  ,WIJ.U  1 1 
of  ifrraraa  Vanda monda  Output  Matrix. 


t taa  papa  AH*  for  flowchart. 


taa  papa  A-44  for  flowchart. 


Figure  A- 21 . Flowchart  for  matrix  element  computation 
steps  for  subroutine  VANINV. 


Specify  Storage  Required  for 
Evaluation  and  Copy  Proper 
Polynomials  for  Numerator  and 
Denominator  into  Working  Arrays. 
— (P/0  VANINVi 


100 y Specify  storage. 
NNS-NON (L) 

NDS-NOD(I) 

For  Numerator  . . 
1 ICOPY  - 1 1 


PNR ( ICOPY) -NR (L , ICOPY) 
PNI ( ICOPY) -N I a , ICOPY) 

-!  ^ Done? 


ICOPY- ICOPY+1 


NNS  > 
wICOPY 


Yes  ...  For  denominator 


ICOPY  - 1 


PDR( ICOPY) -DR (I, ICOPY) 
PDI (ICOPY) =DI(If ICOPY) 

-i  Done? 


IlCOPY— ICOPY-f  1 


TJDS  >. 
ICOPY 


Set  gains. 


RNI-0. 
RDR-DKR(I) 
RDI-DKI (I) 


End  of  Storage  Specification 
and  Polnomial  Copying. 


Figure  A- 22.  Flowchart  of  polynomial  transfer  steps  for 

subroutine  VANINV  (if  no  differentiation  required) 


■*_  ' la  ' i " 
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specify  Storage  Required  for  Differentiation 
and  Copy  Proper  Polynomials  for  ' u*  and  ' v' 
into  Working  Arrays. 

(P/0  VANINV)  — 

11 Gy  Specify  storage. 
NUS*NON  ( L)_ 

NVS=NOD(I). 

NT  =NUS+NVS 
NNS=NT-2 
NDS*  2 *NVS- 1 


For  U ... 


icopy*i 


PUI ( ICOPY) *NI (L, ICOPY) 


Done? 


IC0PY*IC0PY+1 


ICOPY 


Yes  . . . For  v ... 


I COPY* 1 


PVR ( ICOPY) *DR( I , ICOPY) 
PVI  ( ICOPY) *DI ( I , ICOPY) 

, , ,1  Done  ? 


ICOPY* ICOPY+1 


NVS  2. 
ICOPY 


UR-1. 
UI-O. 
VR-DKR(I) 
VI—DKI (I) 


. Set  gains. 


End  Storage  Specification 
and  Polynomial  Copying. 


Figure  A-2  3 . Flowchart  of  polynomial  transfer  steps  for  subroutine 
VANINV  (if  differentiation  required) . 


Find  Derivative 
^ (P/0  VANINVL^ 


Call  FRACDIF(PUR,PUI,PVR,PVI,NUS,NVS,PNR,PNI,PDR, 
PDI,NNS,NDS,NT,UR,UI,VR,VI,RNR,RNI, 
RDR.RDI ,T1R,T1I,T2R,T2I ,T3R,T3I) 


NDIV-NDIV-1 


More  derivatives  needed? 


NDIV  > 0 


Yes  ...  Setup  for 
next  derivative. 

Specify  new 
storage . 


NUS-NNS 
NVS*NDS 
NT  -NUS+NVS 
NNS-NT-2 
NDS=2  *NVS- 1 


End  Fin 
Derivative 


For  new  u ... 


Call  COPYPOL  (PNR,PN I, RNR,RNI ,NUS,PUR,PUI ,UR,UI) 


For  new  v ... 


Cal 1 COPYPOL (PDR, PDI , RDR , RDI , NVS , PVR , PVI , VR, VI ) 


Figure  A-24.  Flowchart  of  differentiation  steps  for 
subroutine  VANINV. 
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V_  Preface  for  Appendix  B 

Two  types  of  programs  are  contained  in  this  appendix.  Part  I 
contains  all  the  subroutines  described  in  Appendix  A.  They  are 
listed  in  alphabetical  order  and  form  a complete  package  of  sub- 
routines for  implementing  Csaki’s  generalized  algorithm  as  described 
in  Section  III  of  the  main  text.  Part  II  contains  programs  and  sub- 
routines needed  to  dimension  the  arrays  for  the  subroutines  of  Part  I . 
They  also  perform  input  and  output  operations,  and  the  tests  that 
are  described  in  Section  IV  of  the  main  text.  Since  the  programs  in 
Part  II  are  not  flowcharted,  a word  description  of  each  is  included 
at  the  beginning  of  Part  II. 
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Descriptions  of  Programs  and  Subroutines  B-30 

Program  READER B-31 

Subroutine  CATEST  B-32 

Subroutine  0UT2 B-35 

Subroutine  UEKTST B-36 

Subroutine  VANF B-37 
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Program  READER  contains  the  dimension  statements  needed  to  test 
subroutine  VANINV' s computational  time  and  accuracy.  After 
reading  in  the  system  eigenvalues,  READER  calls  subroutine 
CATEST  to  perform  the  tests.  In  the  form  shown,  READER  works 
only  with  real,  distinct  eigenvalues.  Some  simple  modifi- 
cations would  enable  it  to  accommodate  complex  and/or  repeated 
eigenvalues . 

Subroutine  CATEST  performs  computation  accuracy  and  time  tests  on 
subroutine  VANINV  and  the  IMSL  subroutine  LINV2F. 

Subroutine  0UT2  prints  out,  in  column  form,  a two  dimensional  array. 

It  is  similar  to  0UT2C  shown  in  Part  I. 

Subroutine  UERTST  disables  the  warning  messages  from  the  IMSL  routines. 
Terminal  errors  CIER  > 1281  are  not  suppressed. 

Subroutine  VANF  forms  the  Vandermonde  matrix  from  a collection  of 
eigenvalues.  In  its  present  form,  it  cannot  handle  complex 
or  repeated  eigenvalues. 
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APPENDIX  C 


SUPPLEMENTARY  MATHEMATICAL  DEVELOPMENTS 


I 


Three  mathematical  developments  are  included  here  to  aid  in  under' 


Augmented  System  State  Transition 


C2 . Determination  of  State  Transition 


Matrices  Using  the  Cayley-Hamilton 


Cl . Sensitivity  Derivatives  as  Part  of  the  Augmented  System 
State  Transition  Matrix 


Using  Eq(8) , the  system  plant  matrix  may  be  augmented  to  find  the 


have  not 


original  n modes,  but  they  are  entirely  independent  of  the  original 


From  Eq  (Cl-11  the  state  transition  matrix  for  the  augmented 


can  be  applied  to  each  section  of  the  partitioned  matrix.  Thus 


which  shows  that  the  sensitivity  derivative  of  the  original  state 


C2.  Determination  of  State  Transition  Matrices 

Using  the  Cayley-Hamilton  Theorem 

This  section  shows  how  the  Cayley-Hamilton  theorem  can  be  applied 
in  order  to  determine  a state  transition  matrix  (STM) . The  use  of  the 
Vandermonde  matrix  and  remainder  polynomial  coefficients  is  described 
for  this  application.  Finally,  the  method  is  extended  to  show  how  it 
can  be  used  to  find  the  partitioned  STM  for  system  sensitivity  calcu- 
lations. Numerical  examples  of  the  concepts  presented  here  are  included 
in  Appendix  D,  Examples  1 and  2. 

Basic  Cayley-Hamilton  Theorem  Application.  The  Cayley-Hamilton 
theorem  states  that  any  square  matrix  satisfies  its  own  characteristic 
equation.  From  it,  the  STM  can  be  shown  to  be  a summation  of  products 
of  the  "remainder  polynomial"  coefficients,  a^Ctl,  and  powers  of  the 
system  matrix.  A,  (Ref  9:V-110  through  V-115b) s 

n 

eAt  - ^ A"3"1  a j Ctl  CC2-U 

j-1 

The  remainder  polynomial  coefficients  can  be  uniquely  determined  from 
n linear  "modal"  equations  CRef  9:V-115bl: 


2-2 


Notice  that  in  Eq  CC2-9)  the  state  transition  matrix  is  expressed  in 
terms  of  products  of:  Cl},  elements  of  the  inverse  transposed  Vander- 
monde matrix,  C21  the  system  modes,  and  (3)  powers  of  the  original 
system  matrix.  The  remainder  polynomial  coefficients,  though  they 
were  used  in  the  theoretical  formulation  of  the  solution,  do  not  even 
appear  in  it  and  thus  it  is  never  necessary  to  actually  calculate 
them.  A very  significant  result  of  this  formulation  is  that  the 

07 


•*  * : -v  v v * 


- _r  \ : * - — • 


original  system  inodes  are  completely  separate  and  visible,  thus 
physical  insight  into  the  system  is  preserved.  This  decoupling  of 
modes  makes  it  possible  to  represent  the  state  transition  matrix  in 
terms  of  the  modes,  each  with  a "component  matrix"  multiplier.  An 
explanation  of  component  matrix  formation  is  given  in  section  C3. 

Extended  Application  of  the  Cayley-Hamilton  Theorem.  The  pro- 
cedure outlined  in  the  preceding  discussion  can  be  very  simply  ex- 
tended to  the  augmented  system  state  transition  matrix.  Eg  (16)  , 
which  is  used  in  the  sensitivity  parameter  calculations  for  the  system 
identification  process.  Eg  (161  is  repeated  here  for  the  reader's 
benefit: 

2n 

*-E 

j«i 

Note  from  Eg  CC2-101  that  the  augmented  system  has  the  same  eigenvalues 
as  the  original  system  since  the  augmented  system  matrix  is  triangular. 
The  only  difference  is  that  the  eigenvalue  multiplicities  of  the  aug- 
mented system  are  double  those  of  the  original  system.  This  means 
that  the  Vandermonde  matrix  for  the  augmented  system  is  2n  * 2n,  there 
are  2n  remainder  polynomial  coefficients,  and  there  are  2n  modes. 
However,  half  of  the  modes  are  simply  derivative  modes  resulting  from 
the  doubled  eigenvalue  multiplicity  which  came  from  augmenting  the 
system.  These  extra  modes  ultimately  have  null  component  matrices 
(Ref  9 : V-1221  and  may  be  thought  of  as  placeholders  or  "pseudo”  modes. 

Thus  it  is  seen  how  the  Cayley-Hamilton  theorem  can  be  applied  to 
determine  the  state  transition  matrices  of  both  the  original  and  aug- 
mented (sensitivity]  systems.  The  next  step  is  to  determine  how  to 


A | 0 
I 

A ' A 
_U)I  _ 


A ! 0 

_ _ J__ 

~ ~~\ 

ACiLl  A 


a^tl 


(C 2-101 
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find  the  "component  matrices"  which  were  revealed  as  part  of  Cayley- 
Hamilton  theorem  application.  That  determination  is  the  subject  of 
the  next  section. 

C3.  Definition  and  Construction  of  Component  Matrices . Several 
references  have  been  made  to  "component  matrices"  in  other  sections  of 
this  thesis.  The  purpose  of  this  section  is  to  show  what  continent 
matrices  are  and  how  they  are  formed. 

A component  matrix  may  be  defined  as  a matrix  which  shows  the 
contribution  of  a given  system  mode  to  the  state  transition  matrix 
(Ref  15:6101.  Thus  it  can  be  seen  from  Eq  (C2-9)  that  each  column 
of  terms  constitutes  a component  matrix,  multiplied  by  its  corres- 
ponding mode. 

Component  matrices  are  commonly  labeled  with  a double  subscript, 

e.g.:  Z^Q.  The  first  subscript  identifies  which  eigenvalue  made  the 

mode  with  which  the  component  matrix  belongs.  The  second  subscript 

can  be  nonzero  only  when  the  system  has  repeated  eigenvalues.  It  then 

identifies  the  number  of  derivatives  which  must  be  taken  of  the 

associated  primary  mode  to  obtain  the  pseudo  mode  with  which  the 

component  matrix  belongs.  (This  derivative  process  is  described  in 

Section  C2.1  Two  examples  should  help  clarify  this.  First  consider 

Z . This  is  the  component  matrix  associated  with  the  primary  mode  of 
10  X.t 

eigenvalue  number  one:  e . Next  consider  Zg^.  This  is  the  component 

matrix  associated  with  one  of  the  pseudo  modes  of  eigenvalue  number 

eight.  Specifically,  it  is  for  the  pseudo  mode  found  by  taking  two 

2 2 ^8^ 

time  derivatives  of  the  primary  mode:  _d e ° t e 

dt2 

Note  that  this  implies  eigenvalue  eight  has  a multiplicity  of  at  least 
three . 

C-9 
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Applying  the  observations  made  and  the  subscript  conventions  just 


described  results  in  a component  matrix  formula: 


n 

j,k  ’Tl/'i.c’'1'1’  -vk+i 


(C3-11 


where  0 < k i ~ 1;  m^  being  the  multiplicity  of  the  j ' th  eigenvalue. 
The  state  transition  matrix  may  then  be  found  from 


M m.-l 

£ £ 


X ,t 


(C3-2) 


in  which  M is  the  number  of  different  eigenvalues.  Some  numerical 
examples  of  the  formation  of  component  matrices  are  given  in 
Appendix  0,  Examples  1 and  2. 


C-10 


iV'l  •**£**?  " 1‘ 


Several  numerical  examples  are  provided  here  to  help  clarify  the 
application  of  concepts  presented  in  the  text.  Examples  1 and  2 work 
with  systems  whose  eigenvalues  were  carefully  selected  for  illustrative 
purposes.  Each  starts  with  a set  of  eigenvalues,  X, , and  a plant 
matrix.  A,  which  contains  a set  of  system  parameters,  6..  For  the 
system  identification  process  outlined  in  the  text,  the  quantities 
sought  are  the  output  sensitivities,  as  defined  by  Eq  Cl 3)  . 

Example  4 illustrates  the  use  of  Csaki's  generalized  algorithm.  This 
is  also  the  process  used  by  subroutine  VANINV. 


- ‘ "■  , ;.  v-  V •• 


Example  1 - Distinct  Eigenvalue  System 
Suppose  n ■ 3,  and 


61--1 


X2  - -3  e2  - -2 


1 


0 

X. 


9. 


0 

0 

X. 


tDl-1) 


The  solution  is  presented  in  four  steps: 

1.  Determination  of  inverse  Vandermonde  matrix  and  (n-1)  powers 
of  the  plant  matrix. 


2.  Calculation  of  the  state  transition  matrix,  e‘ 


At 


3.  Calculation  of  the  parameter  sensitivity  derivatives,  e 


At 


Ci)  ' 


of  the  state  transition  matrix. 


4.  Determination  of  eA  and  eA^  by  means  of  component  matrices 


rather  than  the  solution  of  n linear  equations  as  done  for 
steps  2 and  3. 


Inverse  Vandermonde  Matrix.  For  distinct  eigenvalues,  no  deriva- 
tives need  be  taken  to  find  the  Vandermonde  matrix: 


Powers  of  Plant  Matrix.  Since  n=3,  only  powers  of  A through 


2 

n-1  - 2 need  be  found.  If  A is  multiplied  out,  the  element  in  the 

i'th  row  and  j'th  column  of  the  product  may  be  represented  by 

n=3 


k-1 


(Dl-4) 


where  the  ot's  represent  the  elements  of  the  original  matrix.  For 
this  example: 


— 

— 

" — — 

xx  Q Q 

>• 

o 

o 

X2  0 0 

0X  A2  0 

81  *2  0 

3 

e1(x1+x2)  x2  o 

j>2  0 X3 

o 

M . 

e2cx1+x3i  o x2 

(Dl-5} 


or 


A 


2 


4 0 Tj 

4 0 0 

16  0 0 

-i-3  a 

-1  -3  0 

-19  0 

-2  0 1 | 

-2  0 1 

-10  0 1 

(Dl-6) 


State  Transition  Matrix.  Eq  CC2-8)  for  this  example  is 


eAt  = aiCtlI  + a2(.t)A  + a3(t)A2 


CD1-7) 


The  remainder  polynomial  coefficients,  (tl , may  be  found  by  applying 
Eq  CC2-71  : 


~ 

At 

a^Ctl 

V V V 

11  12  13 

e 1 

Xjt 

a2  Ctl 

m 

v_.  v v.. 

e 

21  22  23 

At 

a3(.t) 

v,,  v, „ 

e 3 

31  32  33 

— — 





(Dl-8) 
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This  results  in 


a1(.tl  * -1/7  e4t  + 1/7  e_3t  + e* 
a2(t)  - 2/21  e4t  - 5/28e"3t  + l/12et 
a3Ctl  - 1/21  e4t  + l/28e”3t  - 1 /Ue* 


(Dl-9) 


Substituting  into  Eq  (D2-71 : 


At  , . 4t  . -3t  t, 

e * (.-1/7  e + 1/7  e + e 1 


0 

1 

0 


+ (2/21  e4t  - 5/28  e“3t  + 1/12 


©jL  -3 

e2  a 


+ (1/21  e4t  + 1/28  e”3t  - 1/12  e*) 


116 


6i  9 

se2  0 


(Dl-10) 


(The  parameter  sensitivities,  9^,  Jure  left  in  for  the  present  so  the 
partials  of  the  state  transition  matrix  w.r.t.  the  parameter  sensi- 
tivities can  be  taken  later.) 

Adding  everything  up  results  in 


At 

e m 


Ce 

e 


4t 


1 4t 


2 4t 


'1  -3t 


2_t 


(Dl-11) 


D-5 


Which,  when  evaluated  at  the  values  of  0^  yields 


C-l/7e4t  + l/7e"3t 


t-2/3e 


+ 2/3e  ) 


e"3t  o (Dl-12) 


Parameter  Sensitivity  Derivatives.  The  parameter  sensitivity 
derivatives,  e*^ , of  the  state  transition  matrix  can  be  found  by 
differentiating  Eq  (pl-11)  w.r.t.  the  parameters,  0^.  For  0s 


Q Q 


At  4t  . “3t 

e^j  ■ Cl/7e  - l/7e  ) 0 0 


0 0 


CD1-131 


And  for  0^: 


0 0 


0 0 


U/3e4t  - l/3efc)  0 0 


CD1-14) 


Use  of  Component  Matrices.  There  exists  one  component  matrix 
for  each  system  mode  or  pseudo  mode,  therefore  it  is  necessary  to  find 
three  component  matrices  for  this  system.  Eq  (C3-2) , which  is  re* 
peated  here  as  Eq  (Dl-15)  may  be  used  to  find  the  state  transition 


matrix. 


M “j"1 


j-1  *-0  j ,k 


k 1 
tZ,  > t*  e 3 


(Dl-15) 


For  this  example,  there  are  M«3  distinct  eigenvalues.  There  are  no 
pseudo  modes,  therefore  k»  0 always.  The  three  component  matrices 


I • 


to  be  found  are  Z^Q,  Z20  Z30'  W*L*-C*1  are  f°un<*  by  using  Eq  IC3-1)  , 

repeated  here  as  (Dl-16) : 


(Dl-16) 


Therefore,  substituting  from  previous  calculations; 


16  0 


9 

0 


0 

0 

1 

0 

0 

0 


Q 

0 

0 


n- 3 


2.-1 


(Dl-1 7a) 


CDl-17bl 


CD1-17  cl 


(Dl-18a) 
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t.t*  *v 
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- . ' 
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x,t  X t X t 

Z,  e + Z e + Z e 
10  20  30 


0 0 


CDl-20b) 


0 0 e + 

la  0 0 


0 0 


0 0 


1 0 


0 0 


^2  Q 1 


0 0 | e 


(Dl-20c) 


<vc 

<h-4t 


1Q  -3t 
70ie 


0 CD1-  20d] 


Comparing  Eq  (JDl-20d)  with.  Eq  CDl-111  shows  that  the  same  answer  comes 
out  either  way.  Note  from  Eq  CDl-20c)  how  clearly  visible  the  modes 
are.  Their  individual  influence  on  the  state  transition  matrix  is  easily 
discernable.  Also,  the  influence  of  the  parameters,  6^,  on  both  the 
modes  and  the  state  transition  matrix  are  easily  seen  from  this 
formulation . 

Note  that  in  actual  practice  the  parameters  are  numbers  rather 
than  variables  as  shown  in  this  example.  They  were  carried  as  variables 
in  the  example  to  show  how  they  influence  the  modes  and  component 
matrices.  Another  important  point  to  note  is  how  direct  the  component 
matrix  approach  turns  out  to  be  in  the  end.  Although  the  development 
may  seem  to  follow  several  independent  paths,  the  application  is  very 
direct.  Once  the  component  matrices  are  found  from  the  inverse  Vander- 
monde matrix,  all  the  effects  of  the  modes  and  parameter  sensitivities 
are  related  to  the  system  by  Eq  (21) . 


v Wi.- 


and 


x2  Q 1 0 0 

1 | 

ei<yy  1 ° ° 

\l  0 J ° 0 

_9L_A2,_!_  0 

o 

o 

M NJ 

o 

0 0 1 X 0 

1 

201(VX2)  ° | ei(VX21  X2 

261  0 1 *1  h 

CD2-7a) 


or 


< 

0 

1 Q 

0 

.j Q 

Q 

0 

0 

1 X1 

Q 

201Q1+X1X2+X21- 

0 

| ejcxj+x^tx*) 

i 

^2 

A3 


4 

9 

0 

-6 


0 | Q 

1 I 0 

Tpr 

0 | 9 


2 

9 

0 

-6 


0 I 0 

±L°_ 

0 | 2 

0 ! -3 


o 

0 

0 

-1 


(D2-7b) 


(J32-8a) 


(D2-8b) 


Notice  from  Eqs  (D2-2,  D2-5,  D2-7)  that  the  parameter  and  its  sensi- 
tivity remain  in  their  respective  positions  for  all  the  powers  of  the 
augmented  matrix.  This  is  an  important  property  because  it  causes  all 
the  modes  to  remain  visible  throughout  the  component  matrix  process  of 
system  identification. 
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State  Transition  Matrix.  For  this  example,  Eq  (C2-8)  is 


eAt  * a^Ct)!  + a2tt}A  + a3(t)A2  + a4Ct)A3  (D2-9) 


Eq  CC2-71  results  in  the  following  system  of  linear  equations,  from 


which  the  remainder  polynomial  coefficients,  a_. (t) , may  be  found: 


ax(tl 
a2  Ct). 
a3Ctl 
a4Ctl 


V V V V 

11  12  13  14 


V21  V22  V23  V24 


V31  V32  V33  V34 
V41  V42  V43  V44 


(D2-10) 


By  applying  Eq  (D2-4) , the  following  coefficients  are  found: 


ax(tl 


a2Ctl 


a3Ctl 


a4Ctl 


21/81e2t  - 18/81te2t  + 60/81e_t  + 36/81te-t 
36/81e2t  - 27/81te2t  - 36/81e“t 


9/8  le 


- 9/81e_t  - 27/81 te_t 


-6/81e2t  + 9/81te2t  + 6/81e_t  + 9/81te_t 


Substituting  into  Eq  (D3-9) : 


eAt  =•  C21/81e2t  - 18/81te2t  + 60/81e_t  + 36/81te_tI 


0)2-11) 


+ C36/81e2t  - 27/81te2t  - 36/81e_t) 


91  -1 


0 

26 


0 

0 
I 

0 | 2 


0 

Q 

0 


i ° ! ei  -i 


+ (9/81e 


2t 


-9/81e-t  - 27/81te“t) 


0 | 0 
1 I 0 


0 

26, 


° I 

0 I 


0 

0 

0 


+ (,-6/81e2t  + 9/81te2t  + 6/81e_t  + 9/81te_t) 


36 


0 

68 


8 0 | Q 

2 -1  I Q 

”oT 


8 


0 

0 

0 


I 2 

i ° I39'  -1 


Summing  up  all  the  terms  yields 


At 
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2t 

e 

a 1 

1 

Q 

0 

2 2t 

0^(27/81e  - 

27/81e”t) 

1 

e_t  j 

0 

Q 

Q 

it 

e 

0 

2t 

-t 

| 

2 2t 

-t  -t 

0 C54/81e"  - 

54/81e  ) 

0 

8^C27/81e  - 

27/81e  1 e 

(D2-13) 

If  Eq  (D2-131  is  evaluated  at  - -3,  the  result  is 
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Applying  these  to  Eq  CD2-151 s 


At  _d 

S(l)  ' d0. 


d0, 


d0. 


d6. 


1 0 

— L 

2 0 


fl 

4 O' 


iJ  i 
8 0 


0 0 


20l  0 

“o  ~ 


20x  Q 
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A comparison  of  Eq  CD2-21)  and  Eq  (D2-13)  shows  that  the  parameter 
sensitivity  of  the  original  state  transition  matrix  is  in  fact  the 
lower  left  block  of  the  partitioned  state  transition  matrix  for  the 
augmented  (sensitivity!  system.  Thus  the  validity  of  Eq  (15)  is 
demonstrated . 

Use  of  Component  Matrices.  Since  the  augmented  system  has  two 
distinct  modes  and  two  pseudo  modes,  a total  of  four  component  matrices 
must  be  found:  2 , Z , Z , Z . The  augmented  system  state  tran- 

J.U  11  fcU  a! 

sition  matrix  and  the  conponent  matrices  may  be  found  using  Eqs  (C 3—  2) 
and  (C3-1) , respectively.  They  are  repeated  here  as  Eqs  (D2-221  and 
(D2-23) . 
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Applying  previous  calculations: 


(D2-22) 


(D2-23) 


Z10  • ]C  Tt.i  A‘~1 


(D2-25a) 


Q 1 

21/81  

0 0 


9!  1 

+ 9/81  — ■ 


— + 36/81 

0 


I - 6/81 


0*  -1  ' 0 0 

0 0 ! 2 0 


• 2 

0 I 91  ‘1 


8 0 | 0 


3B1  -1  I 0 0 

0 0?  8 0 


601  0 i -1 


(D2-25b) 


I 


u 


I 


( 


o 


0 

0 

Q 

0 


Q | 0 

0 I 0 


0 I ° 
0 | 0 


0 

0 

0 

0 


(D2-28C) 


It  can  be  seen  from  the  above  that  if  the  component  matrices  are  cal- 
culated in  order  of  increasing  value  of  the  subscript,  the  calculation 

of  c,  the  column  index  for  v , need  not  be  computed.  It  can  be 

X*  0 c 

simply  initialized  to  one  for  the  first  component  matrix,  and  there- 
after incremented  by  one  for  each  new  component  matrix  calculation. 
Substituting  the  component  matrices  into  Eq  (D2-22J  yields  the 


state  transition  matrix  for  the  augmented  systems 
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Comparing  Eg  (D2-29d)  with.  Eg  (D2-13)  shows  that  the  same  solution  is 
found  by  either  method.  Again  the  influence  of  the  system  modes  and 
parameters  is  clearly  visible  throughout  the  component  matrix  method 
solution  process. 


Example  3 - Use  of  Csaki’s  Generalized  Algorithm 

An  explanation  of  Csaki's  generalized  algorithm  is  found  in 
Section  III  of  the  main  text.  Suppose  for  this  example  that  a fourth 
order  system  has  one  distinct  eigenvalue  and  one  eigenvalue  repeated 
three  times:  X^  * 1 with  multiplicity  k.  ■ 3,  and  X 2 = -1  with 
multiplicity  k ■ 1.  The  Vandermonde  matrix  can  then  be  written 
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which  has  an  inverse 
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The  characteristic  polynomial  is 
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There  are  two  denominator  polynomials — one  for  each  eigenvalue . They 
are 
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For  every  application  of  Csaki’s  algorithm  there  are  n truncated  poly- 
nomials, N£(S) . In  this  example  they  are: 
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The  following  structure  describes  the  inverse  Vandermonde  for  this 
example : 
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Each  element  is  evaluated  as  follows,  using  Eq  (251 , repeated  here  as 
Eq  CD  3-111  : 
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